From 94e753cb8dbfa947b9ecafbae9a703188dc97085 Mon Sep 17 00:00:00 2001 From: Pascal Kaiser Date: Mon, 11 Apr 2022 10:23:11 +0200 Subject: [PATCH 01/38] Add additional type hints --- cara/apps/calculator/model_generator.py | 2 +- cara/apps/calculator/report_generator.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index 595094bc..eab4d069 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -642,7 +642,7 @@ def build_expiration(expiration_definition) -> models._ExpirationBase: return expiration_distribution(tuple(BLO_factors)) -def baseline_raw_form_data(): +def baseline_raw_form_data() -> typing.Dict[str, typing.Union[str, float]]: # Note: This isn't a special "baseline". It can be updated as required. return { 'activity_type': 'office', diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index 258ed6bf..390c81a3 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -96,7 +96,7 @@ def interesting_times(model: models.ExposureModel, approx_n_pts=100) -> typing.L return nice_times -def calculate_report_data(model: models.ExposureModel): +def calculate_report_data(model: models.ExposureModel) -> typing.Dict[str, typing.Any]: times = interesting_times(model) concentrations = [ From 1279ae2549cf1c8b5b1d13ffcfb0f25fafdd0194 Mon Sep 17 00:00:00 2001 From: Pascal Kaiser Date: Mon, 11 Apr 2022 10:24:10 +0200 Subject: [PATCH 02/38] Add additional HTTP POST handler (responds with JSON) - Gets algo input from HTTP POST request body (JSON) - Returns 'report_data' in HTTP response body (JSON) --- cara/apps/calculator/__init__.py | 41 +++++++++++++++++++++++++++++++- 1 file changed, 40 insertions(+), 1 deletion(-) diff --git a/cara/apps/calculator/__init__.py b/cara/apps/calculator/__init__.py index b9572896..5e8390bb 100644 --- a/cara/apps/calculator/__init__.py +++ b/cara/apps/calculator/__init__.py @@ -22,7 +22,7 @@ import tornado.log from . import markdown_tools from . import model_generator -from .report_generator import ReportGenerator +from .report_generator import ReportGenerator, calculate_report_data from .user import AuthenticatedUser, AnonymousUser @@ -129,6 +129,44 @@ class ConcentrationModel(BaseRequestHandler): self.finish(report) +class ConcentrationModelJsonResponse(BaseRequestHandler): + def check_xsrf_cookie(self): + """ + This request handler implements a stateless API that returns report data in JSON format. + Thus, XSRF cookies are disabled by overriding base class implementation of this method with a pass statement. + """ + pass + + async def post(self): + """ + Expects algorithm input in HTTP POST request body in JSON format. + Returns report data (algorithm output) in HTTP POST response body in JSON format. + """ + requested_model_config = json.loads(self.request.body) + if self.settings.get("debug", False): + from pprint import pprint + pprint(requested_model_config) + + try: + form = model_generator.FormData.from_dict(requested_model_config) + except Exception as err: + if self.settings.get("debug", False): + import traceback + print(traceback.format_exc()) + response_json = {'code': 400, 'error': f'Your request was invalid {html.escape(str(err))}'} + self.set_status(400) + await self.finish(json.dumps(response_json)) + return + + executor = loky.get_reusable_executor( + max_workers=self.settings['handler_worker_pool_size'], + timeout=300, + ) + report_data_task = executor.submit(calculate_report_data, form.build_model()) + report_data: dict = await asyncio.wrap_future(report_data_task) + await self.finish(report_data) + + class StaticModel(BaseRequestHandler): async def get(self): form = model_generator.FormData.from_dict(model_generator.baseline_raw_form_data()) @@ -222,6 +260,7 @@ def make_app( (r'/static/(.*)', StaticFileHandler, {'path': static_dir}), (calculator_prefix + r'/?', CalculatorForm), (calculator_prefix + r'/report', ConcentrationModel), + (calculator_prefix + r'/report-json', ConcentrationModelJsonResponse), (calculator_prefix + r'/baseline-model/result', StaticModel), (calculator_prefix + r'/user-guide', ReadmeHandler), (calculator_prefix + r'/static/(.*)', StaticFileHandler, {'path': calculator_static_dir}), From 7cebd355716b29e029f1c668ea4859e64c6a1889 Mon Sep 17 00:00:00 2001 From: Pascal Kaiser Date: Mon, 11 Apr 2022 10:24:23 +0200 Subject: [PATCH 03/38] Add test for added HTTP POST handler --- .../tests/apps/calculator/test_report_json.py | 31 +++++++++++++++++++ 1 file changed, 31 insertions(+) create mode 100644 cara/tests/apps/calculator/test_report_json.py diff --git a/cara/tests/apps/calculator/test_report_json.py b/cara/tests/apps/calculator/test_report_json.py new file mode 100644 index 00000000..a15983df --- /dev/null +++ b/cara/tests/apps/calculator/test_report_json.py @@ -0,0 +1,31 @@ +import json + +import tornado.testing + +import cara.apps.calculator +from cara.apps.calculator import model_generator + +_TIMEOUT = 40. + + +class TestCalculatorJsonResponse(tornado.testing.AsyncHTTPTestCase): + def setUp(self): + super().setUp() + self.http_client.defaults['request_timeout'] = _TIMEOUT + + def get_app(self): + return cara.apps.calculator.make_app() + + @tornado.testing.gen_test(timeout=_TIMEOUT) + def test_json_response(self): + response = yield self.http_client.fetch( + request=self.get_url("/calculator/report-json"), + method="POST", + headers={'content-type': 'application/json'}, + body=json.dumps(model_generator.baseline_raw_form_data()) + ) + self.assertEqual(response.code, 200) + + data = json.loads(response.body) + self.assertIsInstance(data['prob_inf'], float) + self.assertIsInstance(data['expected_new_cases'], float) From ce98ba0026a1a3cc3f54ce0edf4d38b9accab3ee Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 13 Apr 2022 11:16:02 +0200 Subject: [PATCH 04/38] Added beta distribution for interpersonal distances. #254 --- cara/monte_carlo/data.py | 6 +++--- cara/monte_carlo/sampleable.py | 17 +++++++++++++++++ cara/tests/models/test_short_range_model.py | 6 +++--- 3 files changed, 23 insertions(+), 6 deletions(-) diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 2c6268e6..b5e15884 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -5,7 +5,7 @@ import numpy as np from scipy import special as sp import cara.monte_carlo as mc -from cara.monte_carlo.sampleable import LogNormal,LogCustomKernel,CustomKernel,Uniform +from cara.monte_carlo.sampleable import LogNormal,LogCustomKernel,CustomKernel,Uniform, Beta sqrt2pi = np.sqrt(2.*np.pi) @@ -202,5 +202,5 @@ short_range_expiration_distributions = { } -# 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) +# Derived from Fig 8 a) "stand-stand" in https://www.mdpi.com/1660-4601/17/4/1445/htm +short_range_distances = Beta(alpha=1.4342766632654418, beta=27.49916410927064) \ No newline at end of file diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index 27907e49..5a29de81 100644 --- a/cara/monte_carlo/sampleable.py +++ b/cara/monte_carlo/sampleable.py @@ -2,6 +2,7 @@ import typing import numpy as np from sklearn.neighbors import KernelDensity # type: ignore +from scipy.stats import beta import cara.models @@ -130,6 +131,22 @@ class LogCustomKernel(SampleableDistribution): return 10 ** kde_model.sample(n_samples=size)[:, 0] +class Beta(SampleableDistribution): + """ + Defines a Beta distribution parameterized by two positive shape parameters, + denoted by alpha (α) and beta (β), that appear as exponents of the random + variable and control the shape of the distribution. + """ + + def __init__(self, alpha: float, beta: float): + # these are resp. the alpha and beta of the underlying distribution + self.alpha = alpha + self.beta = beta + + def generate_samples(self, size: int) -> float_array_size_n: + return beta.rvs(a = self.alpha, b = self.beta, loc=0.5, scale=7, size=size) + + _VectorisedFloatOrSampleable = typing.Union[ SampleableDistribution, cara.models._VectorisedFloat, ] diff --git a/cara/tests/models/test_short_range_model.py b/cara/tests/models/test_short_range_model.py index a07a656b..4a191334 100644 --- a/cara/tests/models/test_short_range_model.py +++ b/cara/tests/models/test_short_range_model.py @@ -72,9 +72,9 @@ def test_dilution_factor(activity, expected_dilution): @pytest.mark.parametrize( "time, expected_short_range_concentration", [ [8.5, 0.], - [10.5, 15.24806213], - [10.6, 15.24806213], - [11.0, 15.24806213], + [10.5, 7.491771], + [10.6, 7.491771], + [11.0, 7.491771], [12.0, 0.], ] ) From 1d9d4930631dca547f1405b838a07b1dd34de019 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 20 Apr 2022 12:12:47 +0200 Subject: [PATCH 05/38] Changed alpha and beta descriptors --- cara/monte_carlo/data.py | 2 +- cara/monte_carlo/sampleable.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index b5e15884..9a6753b8 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -203,4 +203,4 @@ short_range_expiration_distributions = { # Derived from Fig 8 a) "stand-stand" in https://www.mdpi.com/1660-4601/17/4/1445/htm -short_range_distances = Beta(alpha=1.4342766632654418, beta=27.49916410927064) \ No newline at end of file +short_range_distances = Beta(alpha=0.588715, beta=1.50214) \ No newline at end of file diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index 5a29de81..2aeaf441 100644 --- a/cara/monte_carlo/sampleable.py +++ b/cara/monte_carlo/sampleable.py @@ -144,7 +144,7 @@ class Beta(SampleableDistribution): self.beta = beta def generate_samples(self, size: int) -> float_array_size_n: - return beta.rvs(a = self.alpha, b = self.beta, loc=0.5, scale=7, size=size) + return beta.rvs(a = self.alpha, b = self.beta, loc=0.5, scale=(2 - 0.5), size=size) _VectorisedFloatOrSampleable = typing.Union[ From de5552c236753d9c8e50154e7a3c958b5f5a597a Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Thu, 21 Apr 2022 10:00:59 +0200 Subject: [PATCH 06/38] Updated test values dor beta distribution --- cara/monte_carlo/data.py | 2 +- cara/monte_carlo/sampleable.py | 6 ++++-- cara/tests/models/test_short_range_model.py | 6 +++--- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 9a6753b8..a725b8fa 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -203,4 +203,4 @@ short_range_expiration_distributions = { # Derived from Fig 8 a) "stand-stand" in https://www.mdpi.com/1660-4601/17/4/1445/htm -short_range_distances = Beta(alpha=0.588715, beta=1.50214) \ No newline at end of file +short_range_distances = Beta(alpha=0.588715, beta=1.50214, loc=0.5, scale=1.5) \ No newline at end of file diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index 2aeaf441..d9d4afc7 100644 --- a/cara/monte_carlo/sampleable.py +++ b/cara/monte_carlo/sampleable.py @@ -138,13 +138,15 @@ class Beta(SampleableDistribution): variable and control the shape of the distribution. """ - def __init__(self, alpha: float, beta: float): + def __init__(self, alpha: float, beta: float, loc: float, scale: float): # these are resp. the alpha and beta of the underlying distribution self.alpha = alpha self.beta = beta + self.loc = loc + self.scale = scale def generate_samples(self, size: int) -> float_array_size_n: - return beta.rvs(a = self.alpha, b = self.beta, loc=0.5, scale=(2 - 0.5), size=size) + return beta.rvs(a = self.alpha, b = self.beta, loc=self.loc, scale=self.scale, size=size) _VectorisedFloatOrSampleable = typing.Union[ diff --git a/cara/tests/models/test_short_range_model.py b/cara/tests/models/test_short_range_model.py index 4a191334..481860a8 100644 --- a/cara/tests/models/test_short_range_model.py +++ b/cara/tests/models/test_short_range_model.py @@ -72,9 +72,9 @@ def test_dilution_factor(activity, expected_dilution): @pytest.mark.parametrize( "time, expected_short_range_concentration", [ [8.5, 0.], - [10.5, 7.491771], - [10.6, 7.491771], - [11.0, 7.491771], + [10.5, 8.037883241318065], + [10.6, 8.037883241318065], + [11.0, 8.037883241318065], [12.0, 0.], ] ) From 5b715937d2c98090f7f9c7a46bc0131e3963dc93 Mon Sep 17 00:00:00 2001 From: jdevine Date: Thu, 21 Apr 2022 13:42:55 +0200 Subject: [PATCH 07/38] Implemented Humidity + Temperature half life function from Dabish et al. --- cara/models.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/cara/models.py b/cara/models.py index bfa50139..40e8a422 100644 --- a/cara/models.py +++ b/cara/models.py @@ -439,28 +439,28 @@ class Virus: #: Pre-populated examples of Viruses. types: typing.ClassVar[typing.Dict[str, "Virus"]] - def halflife(self, humidity: _VectorisedFloat) -> _VectorisedFloat: + def halflife(self, humidity: _VectorisedFloat, inside_temp: _VectorisedFloat) -> _VectorisedFloat: # Biological decay (inactivation of the virus in air) - virus # dependent and function of humidity raise NotImplementedError - def decay_constant(self, humidity: _VectorisedFloat) -> _VectorisedFloat: + def decay_constant(self, humidity: _VectorisedFloat, inside_temp: _VectorisedFloat) -> _VectorisedFloat: # Viral inactivation per hour (h^-1) (function of humidity) - return np.log(2) / self.halflife(humidity) + return np.log(2) / self.halflife(humidity, inside_temp) @dataclass(frozen=True) class SARSCoV2(Virus): - def halflife(self, humidity: _VectorisedFloat) -> _VectorisedFloat: + def halflife(self, humidity: _VectorisedFloat, inside_temp: _VectorisedFloat) -> _VectorisedFloat: """ Half-life changes with humidity level. Here is implemented a simple piecewise constant model (for more details see A. Henriques et al, CERN-OPEN-2021-004, DOI: 10.17181/CERN.1GDQ.5Y75) """ - # Taken from Morris et al (https://doi.org/10.7554/eLife.65902) data at T = 22°C and RH = 40 %, - # and from Doremalen et al (https://www.nejm.org/doi/10.1056/NEJMc2004973). - return np.piecewise(humidity, [humidity <= 0.4, humidity > 0.4], [6.43, 1.1]) + # Updated to use the formula from Dabish et al. https://doi.org/10.1080/02786826.2020.1829536 + # with a minimum at hl = 1.1 + return max(1.1, (0.693/(0.16030 + 0.04018((inside_temp-20.615)/10.585)+0.2176((humidity-45.235)/28.665)+0.1))) Virus.types = { @@ -917,9 +917,9 @@ class ConcentrationModel: h = 1.5 # Deposition rate (h^-1) k = (vg * 3600) / h - + #todo: Inside_temp needs to be exposed/added to the room; return ( - k + self.virus.decay_constant(self.room.humidity) + k + self.virus.decay_constant(self.room.humidity, self.room.inside_temp) + self.ventilation.air_exchange(self.room, time) ) From 52300df16b5e3c2a55faa55d6d2e055199db1a60 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Fri, 22 Apr 2022 15:20:20 +0200 Subject: [PATCH 08/38] Changed tests - inside temp is now a property of Room --- cara/apps/calculator/model_generator.py | 5 +- cara/apps/expert.py | 28 +++++----- cara/models.py | 56 ++++++++++--------- .../apps/calculator/test_model_generator.py | 7 +-- cara/tests/conftest.py | 3 +- cara/tests/models/test_concentration_model.py | 2 +- cara/tests/test_known_quantities.py | 28 ++++------ cara/tests/test_monte_carlo.py | 4 +- cara/tests/test_monte_carlo_full_models.py | 9 +-- cara/tests/test_ventilation.py | 10 +--- 10 files changed, 66 insertions(+), 86 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index c4593b49..ec770589 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -244,7 +244,7 @@ class FormData: humidity = 0.3 else: humidity = 0.5 - room = models.Room(volume=volume, humidity=humidity) + room = models.Room(volume=volume, inside_temp=models.PiecewiseConstant((0, 24), (293,)), humidity=humidity) infected_population = self.infected_population() @@ -329,13 +329,11 @@ class FormData: window_interval = always_on outside_temp = self.outside_temp() - inside_temp = models.PiecewiseConstant((0, 24), (293,)) ventilation: models.Ventilation if self.window_type == 'window_sliding': ventilation = models.SlidingWindow( active=window_interval, - inside_temp=inside_temp, outside_temp=outside_temp, window_height=self.window_height, opening_length=self.opening_distance, @@ -344,7 +342,6 @@ class FormData: elif self.window_type == 'window_hinged': ventilation = models.HingedWindow( active=window_interval, - inside_temp=inside_temp, outside_temp=outside_temp, window_height=self.window_height, window_width=self.window_width, diff --git a/cara/apps/expert.py b/cara/apps/expert.py index 84beed27..84690d3c 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -243,19 +243,29 @@ class ModelWidgets(View): def _build_room(self, node): room_volume = widgets.IntSlider(value=node.volume, min=10, max=150) + inside_temp = widgets.IntSlider(value=node.inside_temp.values[0]-273.15, min=15., max=25.) - def on_value_change(change): + def on_volume_change(change): node.volume = change['new'] + def on_insidetemp_change(change): + node.inside_temp.values = (change['new']+273.15,) + # TODO: Link the state back to the widget, not just the other way around. - room_volume.observe(on_value_change, names=['value']) + room_volume.observe(on_volume_change, names=['value']) + inside_temp.observe(on_insidetemp_change, names=['value']) + def on_state_change(): room_volume.value = node.volume + inside_temp.value = node.inside_temp + node.dcs_observe(on_state_change) widget = collapsible( [widget_group( - [[widgets.Label('Room volume (m³)'), room_volume]] + [[widgets.Label('Room volume (m³)'), room_volume], + [widgets.Label('Inside temperature (℃)'), inside_temp], + ] )], title='Specification of workplace', ) @@ -281,7 +291,6 @@ class ModelWidgets(View): def _build_window(self, node) -> WidgetGroup: period = widgets.IntSlider(value=node.active.period, min=0, max=240) interval = widgets.IntSlider(value=node.active.duration, min=0, max=240) - inside_temp = widgets.IntSlider(value=node.inside_temp.values[0]-273.15, min=15., max=25.) def on_period_change(change): node.active.period = change['new'] @@ -289,13 +298,9 @@ class ModelWidgets(View): def on_interval_change(change): node.active.duration = change['new'] - def insidetemp_change(change): - node.inside_temp.values = (change['new']+273.15,) - # TODO: Link the state back to the widget, not just the other way around. period.observe(on_period_change, names=['value']) interval.observe(on_interval_change, names=['value']) - inside_temp.observe(insidetemp_change, names=['value']) outsidetemp_widgets = { 'Fixed': self._build_outsidetemp(node.outside_temp), @@ -327,10 +332,6 @@ class ModelWidgets(View): widgets.Label('Duration of opening (minutes)', layout=auto_width), interval, ), - ( - widgets.Label('Inside temperature (℃)', layout=auto_width), - inside_temp, - ), ( widgets.Label('Outside temperature scheme', layout=auto_width), outsidetemp_w, @@ -485,10 +486,9 @@ class ModelWidgets(View): baseline_model = models.ExposureModel( concentration_model=models.ConcentrationModel( - room=models.Room(volume=75), + room=models.Room(volume=75, inside_temp=models.PiecewiseConstant((0., 24.), (293.15,))), ventilation=models.SlidingWindow( active=models.PeriodicInterval(period=120, duration=15), - inside_temp=models.PiecewiseConstant((0., 24.), (293.15,)), outside_temp=models.PiecewiseConstant((0., 24.), (283.15,)), window_height=1.6, opening_length=0.6, ), diff --git a/cara/models.py b/cara/models.py index 40e8a422..d696067a 100644 --- a/cara/models.py +++ b/cara/models.py @@ -57,15 +57,6 @@ _VectorisedFloat = typing.Union[float, np.ndarray] _VectorisedInt = typing.Union[int, np.ndarray] -@dataclass(frozen=True) -class Room: - #: The total volume of the room - volume: _VectorisedFloat - - #: The humidity in the room (from 0 to 1 - e.g. 0.5 is 50% humidity) - humidity: _VectorisedFloat = 0.5 - - Time_t = typing.TypeVar('Time_t', float, int) BoundaryPair_t = typing.Tuple[Time_t, Time_t] BoundarySequence_t = typing.Union[typing.Tuple[BoundaryPair_t, ...], typing.Tuple] @@ -195,6 +186,18 @@ class PiecewiseConstant: ) +@dataclass(frozen=True) +class Room: + #: The total volume of the room + volume: _VectorisedFloat + + #: The temperature inside the room (Kelvin). + inside_temp: PiecewiseConstant = PiecewiseConstant((0, 24), (293,)) + + #: The humidity in the room (from 0 to 1 - e.g. 0.5 is 50% humidity) + humidity: _VectorisedFloat = 0.5 + + @dataclass(frozen=True) class _VentilationBase: """ @@ -207,7 +210,7 @@ class _VentilationBase: mechanical air exchange through a filter. """ - def transition_times(self) -> typing.Set[float]: + def transition_times(self, room: Room) -> typing.Set[float]: raise NotImplementedError("Subclass must implement") def air_exchange(self, room: Room, time: float) -> _VectorisedFloat: @@ -228,7 +231,7 @@ class Ventilation(_VentilationBase): #: The interval in which the ventilation is active. active: Interval - def transition_times(self) -> typing.Set[float]: + def transition_times(self, room: Room) -> typing.Set[float]: return self.active.transition_times() @@ -243,10 +246,10 @@ class MultipleVentilation(_VentilationBase): """ ventilations: typing.Tuple[_VentilationBase, ...] - def transition_times(self) -> typing.Set[float]: + def transition_times(self, room: Room) -> typing.Set[float]: transitions = set() for ventilation in self.ventilations: - transitions.update(ventilation.transition_times()) + transitions.update(ventilation.transition_times(room)) return transitions def air_exchange(self, room: Room, time: float) -> _VectorisedFloat: @@ -265,9 +268,6 @@ class WindowOpening(Ventilation): #: The interval in which the window is open. active: Interval - #: The temperature inside the room (Kelvin). - inside_temp: PiecewiseConstant - #: The temperature outside of the window (Kelvin). outside_temp: PiecewiseConstant @@ -292,9 +292,10 @@ class WindowOpening(Ventilation): """ raise NotImplementedError("Unknown discharge coefficient") - def transition_times(self) -> typing.Set[float]: - transitions = super().transition_times() - transitions.update(self.inside_temp.transition_times) + def transition_times(self, room: Room) -> typing.Set[float]: + transitions = super().transition_times(room) + print(room.inside_temp.transition_times) + transitions.update(room.inside_temp.transition_times) transitions.update(self.outside_temp.transition_times) return transitions @@ -304,7 +305,7 @@ class WindowOpening(Ventilation): return 0. # Reminder, no dependence on time in the resulting calculation. - inside_temp: _VectorisedFloat = self.inside_temp.value(time) + inside_temp: _VectorisedFloat = room.inside_temp.value(time) outside_temp: _VectorisedFloat = self.outside_temp.value(time) # The inside_temperature is forced to be always at least min_deltaT degree @@ -439,20 +440,20 @@ class Virus: #: Pre-populated examples of Viruses. types: typing.ClassVar[typing.Dict[str, "Virus"]] - def halflife(self, humidity: _VectorisedFloat, inside_temp: _VectorisedFloat) -> _VectorisedFloat: + def halflife(self, humidity: _VectorisedFloat, inside_temp: PiecewiseConstant, time: float) -> _VectorisedFloat: # Biological decay (inactivation of the virus in air) - virus # dependent and function of humidity raise NotImplementedError - def decay_constant(self, humidity: _VectorisedFloat, inside_temp: _VectorisedFloat) -> _VectorisedFloat: + def decay_constant(self, humidity: _VectorisedFloat, inside_temp: PiecewiseConstant, time: float) -> _VectorisedFloat: # Viral inactivation per hour (h^-1) (function of humidity) - return np.log(2) / self.halflife(humidity, inside_temp) + return np.log(2) / self.halflife(humidity, inside_temp, time) @dataclass(frozen=True) class SARSCoV2(Virus): - def halflife(self, humidity: _VectorisedFloat, inside_temp: _VectorisedFloat) -> _VectorisedFloat: + def halflife(self, humidity: _VectorisedFloat, inside_temp: PiecewiseConstant, time: float) -> _VectorisedFloat: """ Half-life changes with humidity level. Here is implemented a simple piecewise constant model (for more details see A. Henriques et al, @@ -460,7 +461,8 @@ class SARSCoV2(Virus): """ # Updated to use the formula from Dabish et al. https://doi.org/10.1080/02786826.2020.1829536 # with a minimum at hl = 1.1 - return max(1.1, (0.693/(0.16030 + 0.04018((inside_temp-20.615)/10.585)+0.2176((humidity-45.235)/28.665)+0.1))) + inside_temp = inside_temp.value(time) + return max(1.1, (0.693/(0.16030 + 0.04018*(((inside_temp-274.15)-20.615)/10.585)+0.02176*((humidity-45.235)/28.665)+0.1))) Virus.types = { @@ -919,7 +921,7 @@ class ConcentrationModel: k = (vg * 3600) / h #todo: Inside_temp needs to be exposed/added to the room; return ( - k + self.virus.decay_constant(self.room.humidity, self.room.inside_temp) + k + self.virus.decay_constant(self.room.humidity, self.room.inside_temp, time) + self.ventilation.air_exchange(self.room, time) ) @@ -950,7 +952,7 @@ class ConcentrationModel: """ state_change_times = {0.} state_change_times.update(self.infected.presence.transition_times()) - state_change_times.update(self.ventilation.transition_times()) + state_change_times.update(self.ventilation.transition_times(self.room)) return sorted(state_change_times) @method_cache diff --git a/cara/tests/apps/calculator/test_model_generator.py b/cara/tests/apps/calculator/test_model_generator.py index 5cc9e409..6f118a8e 100644 --- a/cara/tests/apps/calculator/test_model_generator.py +++ b/cara/tests/apps/calculator/test_model_generator.py @@ -60,7 +60,6 @@ def test_ventilation_slidingwindow(baseline_form: model_generator.FormData): window = models.SlidingWindow( active=models.PeriodicInterval(period=120, duration=10, start=minutes_since_midnight(9 * 60)), - inside_temp=models.PiecewiseConstant((0, 24), (293,)), outside_temp=baseline_window.outside_temp, window_height=1.6, opening_length=0.6, ) @@ -92,7 +91,6 @@ def test_ventilation_hingedwindow(baseline_form: model_generator.FormData): window = models.HingedWindow( active=models.PeriodicInterval(period=120, duration=10, start=minutes_since_midnight(9 * 60)), - inside_temp=models.PiecewiseConstant((0, 24), (293,)), outside_temp=baseline_window.outside_temp, window_height=1.6, window_width=1., opening_length=0.6, ) @@ -106,7 +104,7 @@ def test_ventilation_hingedwindow(baseline_form: model_generator.FormData): def test_ventilation_mechanical(baseline_form: model_generator.FormData): - room = models.Room(75) + room = models.Room(volume=75, inside_temp=models.PiecewiseConstant((0, 24), (293,))) mech = models.HVACMechanical( active=models.PeriodicInterval(period=120, duration=120), q_air_mech=500., @@ -121,7 +119,7 @@ def test_ventilation_mechanical(baseline_form: model_generator.FormData): def test_ventilation_airchanges(baseline_form: model_generator.FormData): - room = models.Room(75) + room = models.Room(75, inside_temp=models.PiecewiseConstant((0, 24), (293,))) airchange = models.AirChange( active=models.PeriodicInterval(period=120, duration=120), air_exch=3., @@ -153,7 +151,6 @@ def test_ventilation_window_hepa(baseline_form: model_generator.FormData): # Now build the equivalent ventilation instance directly, and compare. window = models.SlidingWindow( active=models.PeriodicInterval(period=120, duration=10, start=minutes_since_midnight(9 * 60)), - inside_temp=models.PiecewiseConstant((0, 24), (293,)), outside_temp=baseline_window.outside_temp, window_height=1.6, opening_length=0.6, ) diff --git a/cara/tests/conftest.py b/cara/tests/conftest.py index 3cce3eb3..ef4bff1d 100644 --- a/cara/tests/conftest.py +++ b/cara/tests/conftest.py @@ -8,7 +8,7 @@ import pytest @pytest.fixture def baseline_concentration_model(): model = models.ConcentrationModel( - room=models.Room(volume=75), + room=models.Room(volume=75, inside_temp=models.PiecewiseConstant((0., 24.), (293,))), ventilation=models.AirChange( active=models.SpecificInterval(((0., 24.), )), air_exch=30., @@ -55,7 +55,6 @@ def exposure_model_w_outside_temp_changes(baseline_exposure_model: models.Exposu baseline_exposure_model, { 'concentration_model.ventilation': models.SlidingWindow( active=models.PeriodicInterval(2.2 * 60, 1.8 * 60), - inside_temp=models.PiecewiseConstant((0., 24.), (293,)), outside_temp=cara.data.GenevaTemperatures['Jan'], window_height=1.6, opening_length=0.6, diff --git a/cara/tests/models/test_concentration_model.py b/cara/tests/models/test_concentration_model.py index 432b4d80..3b1d9292 100644 --- a/cara/tests/models/test_concentration_model.py +++ b/cara/tests/models/test_concentration_model.py @@ -26,7 +26,7 @@ def test_concentration_model_vectorisation(override_params): always = models.PeriodicInterval(240, 240) # TODO: This should be a thing on an interval. c_model = models.ConcentrationModel( - models.Room(defaults['volume'], defaults['humidity']), + models.Room(defaults['volume'], models.PiecewiseConstant((0., 24.), (293,)), defaults['humidity']), models.AirChange(always, defaults['air_change']), models.InfectedPopulation( number=1, diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 78c4541f..ed775312 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -19,7 +19,6 @@ def test_no_mask_superspeading_emission_rate(baseline_concentration_model): def baseline_periodic_window(): return models.SlidingWindow( active=models.PeriodicInterval(period=120, duration=15), - inside_temp=models.PiecewiseConstant((0., 24.), (293,)), outside_temp=models.PiecewiseConstant((0., 24.), (283,)), window_height=1.6, opening_length=0.6, ) @@ -27,7 +26,7 @@ def baseline_periodic_window(): @pytest.fixture def baseline_room(): - return models.Room(volume=75) + return models.Room(volume=75, inside_temp=models.PiecewiseConstant((0., 24.), (293,))) @pytest.fixture @@ -131,11 +130,10 @@ def test_periodic_hepa(baseline_periodic_hepa, baseline_room): ], ) def test_multiple_ventilation_HEPA_window(baseline_periodic_hepa, time, expected_value): - room = models.Room(volume=68.) + room = models.Room(volume=68., inside_temp=models.PiecewiseConstant((0., 24.),(293.15,))) tempOutside = models.PiecewiseConstant((0., 1., 2.5),(273.15, 283.15)) - tempInside = models.PiecewiseConstant((0., 24.),(293.15,)) window = models.SlidingWindow(active=models.SpecificInterval([(1 / 60, 24.)]), - inside_temp=tempInside,outside_temp=tempOutside, + outside_temp=tempOutside, window_height=1.,opening_length=0.6) vent = models.MultipleVentilation([window, baseline_periodic_hepa]) npt.assert_allclose(vent.air_exchange(room,time), expected_value, rtol=1e-5) @@ -143,12 +141,12 @@ def test_multiple_ventilation_HEPA_window(baseline_periodic_hepa, time, expected def test_multiple_ventilation_HEPA_window_transitions(baseline_periodic_hepa): tempOutside = models.PiecewiseConstant((0., 1., 2.5),(273.15, 283.15)) - tempInside = models.PiecewiseConstant((0., 24.),(293.15,)) + room = models.Room(68, models.PiecewiseConstant((0., 24.),(293.15,))) window = models.SlidingWindow(active=models.SpecificInterval([(1 / 60, 24.)]), - inside_temp=tempInside,outside_temp=tempOutside, + outside_temp=tempOutside, window_height=1.,opening_length=0.6) vent = models.MultipleVentilation([window, baseline_periodic_hepa]) - assert set(vent.transition_times()) == set([0.0, 1/60, 0.25, 1.0, 2.0, 2.25, + assert set(vent.transition_times(room)) == set([0.0, 1/60, 0.25, 1.0, 2.0, 2.25, 2.5, 4.0, 4.25, 6.0, 6.25, 8.0, 8.25, 10.0, 10.25, 12.0, 12.25, 14.0, 14.25, 16.0, 16.25, 18.0, 18.25, 20.0, 20.25, 22.0, 22.25, 24.]) @@ -188,14 +186,13 @@ def test_multiple_ventilation_HEPA_HVAC_AirChange(volume, expected_value): ) def test_windowopening(time, expected_value): tempOutside = models.PiecewiseConstant((0., 10., 24.),(273.15, 283.15)) - tempInside = models.PiecewiseConstant((0., 24.), (293.15,)) w = models.SlidingWindow( active=models.SpecificInterval([(0., 24.)]), - inside_temp=tempInside,outside_temp=tempOutside, + outside_temp=tempOutside, window_height=1., opening_length=0.6, ) npt.assert_allclose( - w.air_exchange(models.Room(volume=68), time), expected_value, rtol=1e-5 + w.air_exchange(models.Room(volume=68, inside_temp=models.PiecewiseConstant((0., 24.), (293.15, ))), time), expected_value, rtol=1e-5 ) @@ -223,10 +220,9 @@ def build_hourly_dependent_model( outside_temp = temperatures[month] model = models.ConcentrationModel( - room=models.Room(volume=75), + room=models.Room(volume=75, inside_temp=models.PiecewiseConstant((0., 24.), (293, ))), ventilation=models.SlidingWindow( active=models.SpecificInterval(intervals_open), - inside_temp=models.PiecewiseConstant((0., 24.), (293, )), outside_temp=outside_temp, window_height=1.6, opening_length=0.6, ), @@ -246,10 +242,9 @@ def build_hourly_dependent_model( def build_constant_temp_model(outside_temp, intervals_open=((7.5, 8.5),)): model = models.ConcentrationModel( - room=models.Room(volume=75), + room=models.Room(volume=75, inside_temp=models.PiecewiseConstant((0., 24.), (293,))), ventilation=models.SlidingWindow( active=models.SpecificInterval(intervals_open), - inside_temp=models.PiecewiseConstant((0., 24.), (293,)), outside_temp=models.PiecewiseConstant((0., 24.), (outside_temp,)), window_height=1.6, opening_length=0.6, ), @@ -271,7 +266,6 @@ def build_hourly_dependent_model_multipleventilation(month, intervals_open=((7.5 vent = models.MultipleVentilation(( models.SlidingWindow( active=models.SpecificInterval(intervals_open), - inside_temp=models.PiecewiseConstant((0., 24.), (293,)), outside_temp=data.GenevaTemperatures[month], window_height=1.6, opening_length=0.6, ), @@ -281,7 +275,7 @@ def build_hourly_dependent_model_multipleventilation(month, intervals_open=((7.5 ), )) model = models.ConcentrationModel( - room=models.Room(volume=75), + room=models.Room(volume=75, inside_temp=models.PiecewiseConstant((0., 24.), (293,))), ventilation=vent, infected=models.EmittingPopulation( number=1, diff --git a/cara/tests/test_monte_carlo.py b/cara/tests/test_monte_carlo.py index a398e716..b948f4b9 100644 --- a/cara/tests/test_monte_carlo.py +++ b/cara/tests/test_monte_carlo.py @@ -40,10 +40,10 @@ def test_type_annotations(): @pytest.fixture 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)), + room=cara.monte_carlo.Room(volume=cara.monte_carlo.sampleable.Normal(75, 20), + inside_temp=cara.models.PiecewiseConstant((0., 24.), (293,))), ventilation=cara.monte_carlo.SlidingWindow( active=cara.models.PeriodicInterval(period=120, duration=120), - inside_temp=cara.models.PiecewiseConstant((0., 24.), (293,)), outside_temp=cara.models.PiecewiseConstant((0., 24.), (283,)), window_height=1.6, opening_length=0.6, ), diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index b97956b9..32fdd460 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -45,12 +45,11 @@ def shared_office_mc(): Corresponds to the 1st line of Table 4 in https://doi.org/10.1101/2021.10.14.21264988 """ concentration_mc = mc.ConcentrationModel( - room=models.Room(volume=50, humidity=0.5), + room=models.Room(volume=50, inside_temp=models.PiecewiseConstant((0., 24.), (298,)), humidity=0.5), ventilation=models.MultipleVentilation( ventilations=( models.SlidingWindow( active=models.PeriodicInterval(period=120, duration=120), - inside_temp=models.PiecewiseConstant((0., 24.), (298,)), outside_temp=data.GenevaTemperatures['Jun'], window_height=1.6, opening_length=0.2, @@ -88,12 +87,11 @@ def classroom_mc(): Corresponds to the 2nd line of Table 4 in https://doi.org/10.1101/2021.10.14.21264988 """ concentration_mc = mc.ConcentrationModel( - room=models.Room(volume=160, humidity=0.3), + room=models.Room(volume=160, inside_temp=models.PiecewiseConstant((0., 24.), (293,)), humidity=0.3), ventilation=models.MultipleVentilation( ventilations=( models.SlidingWindow( active=models.PeriodicInterval(period=120, duration=120), - inside_temp=models.PiecewiseConstant((0., 24.), (293,)), outside_temp=TorontoTemperatures['Dec'], window_height=1.6, opening_length=0.2, @@ -348,12 +346,11 @@ def test_report_models(mc_model, expected_pi, expected_new_cases, def test_small_shared_office_Geneva(mask_type, month, expected_pi, expected_dose, expected_ER): concentration_mc = mc.ConcentrationModel( - room=models.Room(volume=33, humidity=0.5), + room=models.Room(volume=33, inside_temp=models.PiecewiseConstant((0., 24.), (293,)), humidity=0.5), ventilation=models.MultipleVentilation( ( models.SlidingWindow( active=models.SpecificInterval(((0., 24.),)), - inside_temp=models.PiecewiseConstant((0., 24.), (293,)), outside_temp=data.GenevaTemperatures[month], window_height=1.5, opening_length=0.2, ), diff --git a/cara/tests/test_ventilation.py b/cara/tests/test_ventilation.py index 134acf57..201da2c7 100644 --- a/cara/tests/test_ventilation.py +++ b/cara/tests/test_ventilation.py @@ -11,7 +11,6 @@ from cara import models def baseline_slidingwindow(): return models.SlidingWindow( active=models.SpecificInterval(((0, 4), (5, 9))), - inside_temp=models.PiecewiseConstant((0, 24), (293,)), outside_temp=models.PiecewiseConstant((0, 24), (283,)), window_height=1.6, opening_length=0.6, ) @@ -21,14 +20,13 @@ def baseline_slidingwindow(): def baseline_hingedwindow(): return models.HingedWindow( active=models.SpecificInterval(((0, 4), (5, 9))), - inside_temp=models.PiecewiseConstant((0, 24), (293,)), outside_temp=models.PiecewiseConstant((0, 24), (283,)), window_height=1.6, opening_length=0.6, window_width=1., ) def test_number_of_windows(baseline_slidingwindow): - room = models.Room(75) + room = models.Room(volume=75, inside_temp=models.PiecewiseConstant((0, 24), (293,))) two_windows = dataclasses.replace(baseline_slidingwindow, number_of_windows=2) one_window_exchange = baseline_slidingwindow.air_exchange(room, 1) @@ -63,9 +61,6 @@ def test_hinged_window(baseline_hingedwindow, window_width, {'outside_temp': models.PiecewiseConstant( (0, 2, 3), (np.array([20, 30, 28]), np.array([25, 30, 27])) )}, - {'inside_temp': models.PiecewiseConstant( - (0, 20), (np.array([20, 30, 25]), ) - )}, ] ) def test_hinged_window_vectorisation(override_params): @@ -73,11 +68,10 @@ def test_hinged_window_vectorisation(override_params): 'window_height': 0.15, 'window_width': 0.15, 'opening_length': 0.15, - 'inside_temp': models.PiecewiseConstant((0, 2, 3), (20, 25)), 'outside_temp': models.PiecewiseConstant((0, 2, 3), (10, 15)), } defaults.update(override_params) - room = models.Room(volume=75) + room = models.Room(volume=75, inside_temp=models.PiecewiseConstant((0, 2, 3), (20, 25))) t = 0.5 window = models.HingedWindow(models.PeriodicInterval(60, 30), **defaults) if {'window_height', 'opening_length', 'window_width'}.intersection(override_params): From f99dcc5e0cefa7261ec454b9c0a84aa160c28639 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 22 Apr 2022 16:50:17 +0200 Subject: [PATCH 09/38] Fix wrong emission rate in test_monte_carlo_full_models.py::test_small_shared_office_Geneva; decrease tolerance and number of samples subsequently --- cara/tests/test_monte_carlo_full_models.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index b97956b9..caed9c52 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 = 600_000 -TOLERANCE = 0.06 +SAMPLE_SIZE = 500_000 +TOLERANCE = 0.05 # Load the weather data (temperature in kelvin) for Toronto. toronto_coordinates = (43.667, 79.400) @@ -342,7 +342,7 @@ def test_report_models(mc_model, expected_pi, expected_new_cases, ["No mask", "Jul", 9.52, 9.920, 809], ["Type I", "Jul", 1.7, 0.913, 149], ["FFP2", "Jul", 0.51, 0.239, 149], - ["Type I", "Feb", 0.57, 0.272, 162], + ["Type I", "Feb", 0.57, 0.272, 149], ], ) def test_small_shared_office_Geneva(mask_type, month, expected_pi, From ef84f28ab6d298b52d6dc8f0940f89f716ac8031 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 22 Apr 2022 17:13:15 +0200 Subject: [PATCH 10/38] Adding a method to extract the right interaction time boundaries between two given times in the short-range model (cosmetics) --- cara/models.py | 52 ++++++++++++++++++++++++++++++++------------------ 1 file changed, 33 insertions(+), 19 deletions(-) diff --git a/cara/models.py b/cara/models.py index bfa50139..8097fd71 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1168,19 +1168,44 @@ class ShortRangeModel: # 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): + @method_cache + def extract_between_bounds(self, time1: float, time2: float) -> typing.Tuple[float,float]: + """ + Extract the bounds of the interval resulting from the + intersection of [time1, time2] and the presence interval. + If [time1, time2] has nothing common to the presence interval, + we return (0, 0). + Raise an error if time1 and time2 are not in ascending order. + """ + if time1>time2: + raise ValueError("time1 must be less or equal to time2") + + start, stop = self.presence.boundaries()[0] + if (stop < time1) or (start > time2): + return (0, 0) + elif start <= time1 and time2<= stop: + return time1, time2 + elif start <= time1 and stop < time2: + return time1, stop + elif time1 < start and time2 <= stop: + return start, time2 + elif time1 <= start and stop < time2: + return start, stop + + 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] + start, stop = self.extract_between_bounds(time1, time2) jet_origin = self.expiration.jet_origin_concentration() dilution = self.dilution_factor() total_normed_concentration_diluted = ( - concentration_model.integrated_concentration(start_bound, - stop_bound)/dilution/ + concentration_model.integrated_concentration(start, stop + )/dilution/ concentration_model.virus.viral_load_in_sputum ) total_normed_concentration_interpolated = np.interp( @@ -1188,7 +1213,7 @@ class ShortRangeModel: concentration_model.infected.particle.diameter, total_normed_concentration_diluted ) - return (jet_origin/dilution * (stop_bound - start_bound) + return (jet_origin/dilution * (stop - start) ) - total_normed_concentration_interpolated @@ -1297,20 +1322,9 @@ class ExposureModel: """ 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) + start, stop = interaction.extract_between_bounds(time1, time2) + short_range_exposure = interaction.normed_exposure_between_bounds( + self.concentration_model, start, stop) fdep = interaction.expiration.particle.fraction_deposited(evaporation_factor=1.0) diameter = interaction.expiration.particle.diameter From b11dd90e928f316701a41113b3555c353b25a5c4 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 22 Apr 2022 22:21:15 +0200 Subject: [PATCH 11/38] Adding full algo tests with distributions, in particular on the concentration/dose/probability of infection --- cara/tests/test_full_algorithm.py | 247 ++++++++++++++++++++++++------ 1 file changed, 197 insertions(+), 50 deletions(-) diff --git a/cara/tests/test_full_algorithm.py b/cara/tests/test_full_algorithm.py index fb6ade81..6726f3a4 100644 --- a/cara/tests/test_full_algorithm.py +++ b/cara/tests/test_full_algorithm.py @@ -244,7 +244,7 @@ class SimpleShortRangeModel: return dilution - @method_cache +# @method_cache def jet_concentration(self,conc_model: SimpleConcentrationModel) -> _VectorisedFloat: """ virion concentration at the origin of the jet (close to @@ -267,13 +267,13 @@ class SimpleShortRangeModel: def concentration(self, conc_model: SimpleConcentrationModel, time: float) -> _VectorisedFloat: """ compute the short-range part of the concentration, and add it - to the background concentration + to the long-range concentration """ if self.interaction_interval.triggered(time): - background_concentration = conc_model.concentration(time) + lr_concentration = conc_model.concentration(time) S = self.dilution_factor() return (self.jet_concentration(conc_model) - - background_concentration) / S + - lr_concentration) / S else: return 0. @@ -352,8 +352,17 @@ class SimpleExposureModel(SimpleConcentrationModel): epsabs=0.,limit=500)[0] * self.viral_load * self.breathing_rate) + def total_concentration(self, t: float): + """ + total concentration at time t + """ + res = self.concentration(t) + for sr_mod in self.sr_models: + res += sr_mod.concentration(self,t) + return res + @method_cache - def integrated_background_concentration(self,t1: float,t2: float, + def integrated_longrange_concentration(self,t1: float,t2: float, evaporation: float) -> _VectorisedFloat: """ background (long-range) concentration integrated from t1 to t2 @@ -417,7 +426,7 @@ class SimpleExposureModel(SimpleConcentrationModel): epsabs=0.,limit=500)[0] * self.viral_load * 1e-6 * (t2-t1) ) result += sr_model.breathing_rate * ( - res-self.integrated_background_concentration(t1,t2,evaporation) + res-self.integrated_longrange_concentration(t1,t2,evaporation) )/sr_model.dilution_factor() return result @@ -429,7 +438,7 @@ class SimpleExposureModel(SimpleConcentrationModel): """ result = 0. for t1,t2 in self.infected_presence.boundaries(): - result += (self.integrated_background_concentration(t1,t2,self.evaporation) + result += (self.integrated_longrange_concentration(t1,t2,self.evaporation) * self.breathing_rate) result += self.integrated_shortrange_concentration() @@ -468,6 +477,25 @@ def c_model() -> mc.ConcentrationModel: ).build_model(SAMPLE_SIZE) +@pytest.fixture +def c_model_distr() -> 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=virus_distributions['SARS_CoV_2_DELTA'], + mask=models.Mask.types['No mask'], + activity=activity_distributions['Seated'], + expiration=expiration_distributions['Breathing'], + host_immunity=0., + ), + evaporation_factor=0.3, + ).build_model(SAMPLE_SIZE) + + @pytest.fixture def sr_models() -> typing.Tuple[mc.ShortRangeModel, ...]: return ( @@ -516,10 +544,107 @@ def simple_sr_models() -> typing.Tuple[SimpleShortRangeModel, ...]: ) +@pytest.fixture +def expo_sr_model(c_model,sr_models) -> mc.ExposureModel: + return mc.ExposureModel( + concentration_model=c_model, + short_range=sr_models, + exposed=mc.Population( + number=1, + presence=presence, + mask=models.Mask.types['No mask'], + activity=models.Activity.types['Seated'], + host_immunity=0., + ), + ).build_model(SAMPLE_SIZE) + + +@pytest.fixture +def simple_expo_sr_model(simple_sr_models) -> SimpleExposureModel: + return SimpleExposureModel( + 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'], + finf = models.Virus.types['SARS_CoV_2_DELTA'].viable_to_RNA_ratio, + HI = 0., + ID50 = models.Virus.types['SARS_CoV_2_DELTA'].infectious_dose, + transmissibility = models.Virus.types['SARS_CoV_2_DELTA'].transmissibility_factor, + sr_models = simple_sr_models, + ) + + +@pytest.fixture +def expo_sr_model_distr(c_model_distr) -> mc.ExposureModel: + return mc.ExposureModel( + concentration_model=c_model_distr, + short_range=( + mc.ShortRangeModel( + expiration = short_range_expiration_distributions['Breathing'], + activity = activity_distributions['Seated'], + presence = interaction_intervals[0], + distance = short_range_distances, + ).build_model(SAMPLE_SIZE), + mc.ShortRangeModel( + expiration = short_range_expiration_distributions['Speaking'], + activity = activity_distributions['Seated'], + presence = interaction_intervals[1], + distance = short_range_distances, + ).build_model(SAMPLE_SIZE), + ), + exposed=mc.Population( + number=1, + presence=presence, + mask=models.Mask.types['No mask'], + activity=models.Activity.types['Seated'], + host_immunity=0., + ), + ).build_model(SAMPLE_SIZE) + + +@pytest.fixture +def simple_expo_sr_model_distr(c_model_distr) -> SimpleExposureModel: + return SimpleExposureModel( + infected_presence = presence, + viral_load = virus_distributions['SARS_CoV_2_DELTA' + ].build_model(SAMPLE_SIZE).viral_load_in_sputum, + breathing_rate = activity_distributions['Seated'].build_model( + SAMPLE_SIZE).exhalation_rate, + room_volume = 50., + lambda_ventilation= 1., + BLO_factors = expiration_BLO_factors['Breathing'], + finf = virus_distributions['SARS_CoV_2_DELTA' + ].build_model(SAMPLE_SIZE).viable_to_RNA_ratio, + HI = 0., + ID50 = virus_distributions['SARS_CoV_2_DELTA' + ].build_model(SAMPLE_SIZE).infectious_dose, + transmissibility = virus_distributions['SARS_CoV_2_DELTA' + ].transmissibility_factor, + sr_models = ( + SimpleShortRangeModel( + interaction_interval = interaction_intervals[0], + distance = short_range_distances.generate_samples(SAMPLE_SIZE), + breathing_rate = activity_distributions['Seated'].build_model( + SAMPLE_SIZE).exhalation_rate, + BLO_factors = expiration_BLO_factors['Breathing'], + ), + SimpleShortRangeModel( + interaction_interval = interaction_intervals[1], + distance = short_range_distances.generate_samples(SAMPLE_SIZE), + breathing_rate = activity_distributions['Seated'].build_model( + SAMPLE_SIZE).exhalation_rate, + BLO_factors = expiration_BLO_factors['Speaking'], + ) + ), + ) + + @pytest.mark.parametrize( "time", np.linspace(8.5,17.5,12), ) -def test_background_concentration(time,c_model,simple_c_model): +def test_longrange_concentration(time,c_model,simple_c_model): npt.assert_allclose( c_model.concentration(time).mean(), simple_c_model.concentration(time), rtol=TOLERANCE @@ -542,7 +667,7 @@ def test_shortrange_concentration(time,c_model,simple_c_model, ) -def test_background_exposure(c_model): +def test_longrange_exposure(c_model): simple_expo_model = SimpleExposureModel( infected_presence = presence, viral_load = models.Virus.types['SARS_CoV_2_DELTA'].viral_load_in_sputum, @@ -577,7 +702,27 @@ def test_background_exposure(c_model): ) -def test_background_exposure_with_distributions(): +@pytest.mark.parametrize( + "time", [11., 12.5, 17.] +) +def test_longrange_concentration_with_distributions(c_model_distr,time): + simple_expo_model = SimpleConcentrationModel( + infected_presence = presence, + viral_load = virus_distributions['SARS_CoV_2_DELTA' + ].build_model(SAMPLE_SIZE).viral_load_in_sputum, + breathing_rate = activity_distributions['Seated'].build_model( + SAMPLE_SIZE).exhalation_rate, + room_volume = 50., + lambda_ventilation= 1., + BLO_factors = expiration_BLO_factors['Breathing'], + ) + npt.assert_allclose( + c_model_distr.concentration(time).mean(), + simple_expo_model.concentration(time).mean(), rtol=TOLERANCE + ) + + +def test_longrange_exposure_with_distributions(c_model_distr): simple_expo_model = SimpleExposureModel( infected_presence = presence, viral_load = virus_distributions['SARS_CoV_2_DELTA' @@ -597,21 +742,7 @@ def test_background_exposure_with_distributions(): sr_models = (), ) expo_model = mc.ExposureModel( - concentration_model=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=virus_distributions['SARS_CoV_2_DELTA'], - mask=models.Mask.types['No mask'], - activity=activity_distributions['Seated'], - expiration=expiration_distributions['Breathing'], - host_immunity=0., - ), - evaporation_factor=0.3, - ), + concentration_model=c_model_distr, short_range=(), exposed=mc.Population( number=1, @@ -631,31 +762,21 @@ def test_background_exposure_with_distributions(): ) -def test_exposure_with_shortrange(c_model,sr_models,simple_sr_models): - simple_expo_sr_model = SimpleExposureModel( - 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'], - finf = models.Virus.types['SARS_CoV_2_DELTA'].viable_to_RNA_ratio, - HI = 0., - ID50 = models.Virus.types['SARS_CoV_2_DELTA'].infectious_dose, - transmissibility = models.Virus.types['SARS_CoV_2_DELTA'].transmissibility_factor, - sr_models = simple_sr_models, - ) - expo_sr_model = mc.ExposureModel( - concentration_model=c_model, - short_range=sr_models, - exposed=mc.Population( - number=1, - presence=presence, - mask=models.Mask.types['No mask'], - activity=models.Activity.types['Seated'], - host_immunity=0., - ), - ).build_model(SAMPLE_SIZE) +# tests on the concentration with short-range should be skipped until +# one finds a way to avoid the large variability of the concentration +# with short-range 'Speaking' or 'Shouting' interactions +@pytest.mark.skip +@pytest.mark.parametrize( + "time", [10.75, 14.75, 16.] +) +def test_concentration_with_shortrange(expo_sr_model,simple_expo_sr_model,time): + npt.assert_allclose( + expo_sr_model.concentration(time).mean(), + simple_expo_sr_model.total_concentration(time).mean(), rtol=TOLERANCE + ) + + +def test_exposure_with_shortrange(expo_sr_model,simple_expo_sr_model): npt.assert_allclose( expo_sr_model.deposited_exposure().mean(), simple_expo_sr_model.dose().mean(), rtol=TOLERANCE @@ -665,3 +786,29 @@ def test_exposure_with_shortrange(c_model,sr_models,simple_sr_models): simple_expo_sr_model.probability_infection().mean(), rtol=TOLERANCE ) + +@pytest.mark.skip +@pytest.mark.parametrize( + "time", [10.75, 14.75, 16.] +) +def test_concentration_with_shortrange_and_distributions( + expo_sr_model_distr,simple_expo_sr_model_distr,time): + npt.assert_allclose( + expo_sr_model_distr.concentration(time).mean(), + simple_expo_sr_model_distr.total_concentration(time).mean(), + rtol=TOLERANCE + ) + + +def test_exposure_with_shortrange_and_distributions(expo_sr_model_distr, + simple_expo_sr_model_distr): + npt.assert_allclose( + expo_sr_model_distr.deposited_exposure().mean(), + simple_expo_sr_model_distr.dose().mean(), rtol=0.05 + ) + npt.assert_allclose( + expo_sr_model_distr.infection_probability().mean(), + simple_expo_sr_model_distr.probability_infection().mean(), + rtol=0.03 + ) + From 9a4790d407ead266ead34f6d61f1cecd6d12d9b1 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 22 Apr 2022 22:23:49 +0200 Subject: [PATCH 12/38] Decoupling jet exposure and the part subtracted from the background concentration, in the short-range model, in order to treat correctly the diameter integrals. This correct the probability of infection. --- cara/models.py | 70 +++++++++++++++++++++++++++++++++++--------------- 1 file changed, 49 insertions(+), 21 deletions(-) diff --git a/cara/models.py b/cara/models.py index 8097fd71..cbb547c9 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1192,29 +1192,46 @@ class ShortRangeModel: elif time1 <= start and stop < time2: return start, stop - def normed_exposure_between_bounds(self, concentration_model: ConcentrationModel, - time1: float, time2: float): + def _normed_jet_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. + Get the part of the integrated short-range concentration of + viruses in the air, between the times start and stop, coming + from the jet concentration, normalized by the viral load, and + without dilution. """ start, stop = self.extract_between_bounds(time1, time2) - jet_origin = self.expiration.jet_origin_concentration() - dilution = self.dilution_factor() + return jet_origin * (stop - start) - total_normed_concentration_diluted = ( - concentration_model.integrated_concentration(start, stop - )/dilution/ - concentration_model.virus.viral_load_in_sputum + def _normed_interpolated_longrange_exposure_between_bounds( + self, concentration_model: ConcentrationModel, + time1: float, time2: float): + """ + Get the part of the integrated short-range concentration due + to the background concentration, normalized by the viral load + and the breathing rate, and without dilution. + One needs to interpolate the integrated long-range concentration + for the particle diameters defined here. + TODO: make sure any potential extrapolation has a + negligible effect. + """ + start, stop = self.extract_between_bounds(time1, time2) + if stop<=start: + return 0. + + normed_int_concentration = ( + concentration_model.integrated_concentration(start, stop) + /concentration_model.virus.viral_load_in_sputum + /concentration_model.infected.activity.exhalation_rate ) - total_normed_concentration_interpolated = np.interp( + normed_int_concentration_interpolated = np.interp( self.expiration.particle.diameter, concentration_model.infected.particle.diameter, - total_normed_concentration_diluted + normed_int_concentration ) - return (jet_origin/dilution * (stop - start) - ) - total_normed_concentration_interpolated + return normed_int_concentration_interpolated @dataclass(frozen=True) @@ -1292,7 +1309,7 @@ class ExposureModel: # 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). + # lead to wrong results for the probability of infection). dep_exposure_integrated = np.array(self._long_range_normed_exposure_between_bounds(time1, time2) * aerosols * fdep).mean() @@ -1323,24 +1340,35 @@ class ExposureModel: deposited_exposure = 0. for interaction in self.short_range: start, stop = interaction.extract_between_bounds(time1, time2) - short_range_exposure = interaction.normed_exposure_between_bounds( + short_range_jet_exposure = interaction._normed_jet_exposure_between_bounds( self.concentration_model, start, stop) + short_range_lr_exposure = interaction._normed_interpolated_longrange_exposure_between_bounds( + self.concentration_model, start, stop) + dilution = interaction.dilution_factor() 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. + # 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() + # lead to wrong results for the probability of infection). + deposited_exposure += (np.array(short_range_jet_exposure + * fdep).mean()/dilution + - np.array(short_range_lr_exposure * fdep).mean() + * self.concentration_model.infected.activity.exhalation_rate + /dilution) 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 + deposited_exposure += (short_range_jet_exposure + * fdep/dilution + - short_range_lr_exposure * fdep + * self.concentration_model.infected.activity.exhalation_rate + /dilution) # multiply by the (diameter-independent) inhalation rate deposited_exposure *= interaction.activity.inhalation_rate From a86d9a838a766cb670fdd8e3156b0fe34ddaa445 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Fri, 22 Apr 2022 22:25:50 +0200 Subject: [PATCH 13/38] In test_short_range_model, removing test on short-range normed exposure (method does not exist anymore), and adding tests on extract_between_bounds method --- cara/tests/models/test_short_range_model.py | 48 +++++++++++++-------- 1 file changed, 29 insertions(+), 19 deletions(-) diff --git a/cara/tests/models/test_short_range_model.py b/cara/tests/models/test_short_range_model.py index a07a656b..b1dee050 100644 --- a/cara/tests/models/test_short_range_model.py +++ b/cara/tests/models/test_short_range_model.py @@ -45,7 +45,8 @@ def test_short_range_model_ndarray(concentration_model, short_range_model): 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._normed_jet_exposure_between_bounds(concentration_model, 10.75, 10.85), np.ndarray) + assert isinstance(model._normed_interpolated_longrange_exposure_between_bounds(concentration_model, 10.75, 10.85), np.ndarray) assert isinstance(model.short_range_concentration(concentration_model, 14.0), float) @@ -69,6 +70,32 @@ def test_dilution_factor(activity, expected_dilution): ) +def test_extract_between_bounds_raise_on_wrong_order(short_range_model): + model = short_range_model.build_model(1) + with pytest.raises(ValueError, match='time1 must be less or equal to time2'): + model.extract_between_bounds(11.,10.) + + +@pytest.mark.parametrize( + "time1, time2, expected_start, expected_stop", [ + [10., 12., 10.5, 11.], + [10., 10.7, 10.5, 10.7], + [10., 10.45, 0., 0.], + [11.01, 11.5, 0., 0.], + [10.8, 10.9, 10.8, 10.9], + [10.8, 11.5, 10.8, 11.], + [10.5, 11., 10.5, 11.], + ] +) +def test_extract_between_bounds(short_range_model, time1, time2, + expected_start, expected_stop): + model = short_range_model.build_model(1) + np.testing.assert_equal( + model.extract_between_bounds(time1, time2), + (expected_start, expected_stop), + ) + + @pytest.mark.parametrize( "time, expected_short_range_concentration", [ [8.5, 0.], @@ -83,23 +110,6 @@ def test_short_range_concentration(time, expected_short_range_concentration, con 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 + expected_short_range_concentration, rtol=0.02 ) -@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 - ) From 6d9210939d2fee299cc78b3d5dcd77db68f0d0ff Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Sat, 23 Apr 2022 09:45:03 +0200 Subject: [PATCH 14/38] Updating version number --- cara/apps/calculator/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/apps/calculator/__init__.py b/cara/apps/calculator/__init__.py index e7c1d5d6..8d4e5e45 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.1.1" +__version__ = "4.1.2" class BaseRequestHandler(RequestHandler): From a85a05316f6b58f50660a47df823e677bbf0c13e Mon Sep 17 00:00:00 2001 From: jdevine Date: Mon, 25 Apr 2022 11:36:49 +0200 Subject: [PATCH 15/38] Adjusted tests, modified expected outputs for new temperature and humidity dependent halflife method --- cara/models.py | 2 +- cara/tests/models/test_exposure_model.py | 12 ++++++------ cara/tests/test_full_algorithm.py | 10 ++++++++-- cara/tests/test_known_quantities.py | 12 ++++++------ cara/tests/test_monte_carlo_full_models.py | 16 ++++++++-------- 5 files changed, 29 insertions(+), 23 deletions(-) diff --git a/cara/models.py b/cara/models.py index d696067a..1adb810e 100644 --- a/cara/models.py +++ b/cara/models.py @@ -462,7 +462,7 @@ class SARSCoV2(Virus): # Updated to use the formula from Dabish et al. https://doi.org/10.1080/02786826.2020.1829536 # with a minimum at hl = 1.1 inside_temp = inside_temp.value(time) - return max(1.1, (0.693/(0.16030 + 0.04018*(((inside_temp-274.15)-20.615)/10.585)+0.02176*((humidity-45.235)/28.665)+0.1))) + return max(1.1, (0.693/(0.16030 + 0.04018*(((inside_temp-273.15)-20.615)/10.585)+0.02176*((humidity-45.235)/28.665)+0.1))) Virus.types = { diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 44979024..abc476d0 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -179,12 +179,12 @@ def sr_model(): @pytest.mark.parametrize( ["exposed_time_interval", "expected_deposited_exposure"], [ - [(0., 1.), 45.6008710], - [(1., 1.01), 0.5280401], - [(1.01, 1.02), 0.51314096385], - [(12., 12.01), 0.016255813185], - [(12., 24.), 645.63619275], - [(0., 24.), 700.7322474], + [(0., 1.), 48.19316], + [(1., 1.01), 0.566368], + [(1.01, 1.02), 0.551401], + [(12., 12.01), 0.016278], + [(12., 24.), 691.21381], + [(0., 24.), 750.258043], ] ) def test_exposure_model_integral_accuracy(exposed_time_interval, diff --git a/cara/tests/test_full_algorithm.py b/cara/tests/test_full_algorithm.py index fb6ade81..71f10f05 100644 --- a/cara/tests/test_full_algorithm.py +++ b/cara/tests/test_full_algorithm.py @@ -47,6 +47,9 @@ class SimpleConcentrationModel: #: room volume (m^3) room_volume: _VectorisedFloat + #: The temperature inside the room (Kelvin). + #inside_temp: float = 293.15 + #: ventilation rate (air changes per hour) - including HEPA lambda_ventilation: _VectorisedFloat @@ -84,8 +87,11 @@ class SimpleConcentrationModel: """ 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) ) + + return (self.lambda_ventilation + + ln2/(max(1.1, (0.693 / (0.16030 + 0.04018 * (((22) - 20.615) / 10.585) + 0.02176 * ( + (self.humidity - 45.235) / 28.665) + 0.1))))) + #6.43 if self.humidity<=0.4 else 1.1) ) @method_cache def deposition_removal_coefficient(self) -> float: diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index ed775312..4eb94d0b 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -43,7 +43,7 @@ def test_concentrations(baseline_concentration_model): 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], + [0.000000e+00, 2.108144e+01, 1.004741e-12, 2.108144e+01, 4.788592e-26], rtol=1e-6 ) @@ -94,7 +94,7 @@ def test_r0(baseline_exposure_model): # expected r0 was computed with a trapezoidal integration, using # a mesh of 100'000 pts per exposed presence interval. r0 = baseline_exposure_model.reproduction_number() - npt.assert_allclose(r0, 776.941990) + npt.assert_allclose(r0, 781.293793) def test_periodic_window(baseline_periodic_window, baseline_room): @@ -381,8 +381,8 @@ def build_exposure_model(concentration_model, short_range_model): @pytest.mark.parametrize( "month, expected_deposited_exposure", [ - ['Jan', 377.440565819], - ['Jun', 1721.03336729], + ['Jan', 392.994454], + ['Jun', 2127.82386], ], ) def test_exposure_hourly_dep(month,expected_deposited_exposure, baseline_sr_model): @@ -402,8 +402,8 @@ def test_exposure_hourly_dep(month,expected_deposited_exposure, baseline_sr_mode @pytest.mark.parametrize( "month, expected_deposited_exposure", [ - ['Jan', 383.339206111], - ['Jun', 1799.17597184], + ['Jan', 394.000261], + ['Jun', 2239.777906], ], ) def test_exposure_hourly_dep_refined(month,expected_deposited_exposure, baseline_sr_model): diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index 32fdd460..58f0a57d 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -310,13 +310,13 @@ def waiting_room_mc(): @pytest.mark.parametrize( "mc_model, expected_pi, expected_new_cases, expected_dose, expected_ER", [ - ["shared_office_mc", 6.03, 0.18, 3.198, 809], + ["shared_office_mc", 6.524158, 0.195725, 3.413983, 809], ["classroom_mc", 9.5, 1.85, 9.478, 5624], ["ski_cabin_mc", 16.0, 0.5, 17.315, 7966], - ["skagit_chorale_mc",65.7, 40.0, 102.213, 190422], - ["bus_ride_mc", 12.0, 8.0, 7.65, 5419], + ["skagit_chorale_mc",69.84901, 40.0, 121.265911, 190422], + ["bus_ride_mc", 13.311721, 8.918853, 8.6662, 5419], ["gym_mc", 0.45, 0.13, 0.208, 1145], - ["waiting_room_mc", 1.59, 0.22, 0.821, 737], + ["waiting_room_mc", 1.854373, 0.259612, 0.99173, 737], ] ) def test_report_models(mc_model, expected_pi, expected_new_cases, @@ -337,10 +337,10 @@ def test_report_models(mc_model, expected_pi, expected_new_cases, @pytest.mark.parametrize( "mask_type, month, expected_pi, expected_dose, expected_ER", [ - ["No mask", "Jul", 9.52, 9.920, 809], - ["Type I", "Jul", 1.7, 0.913, 149], - ["FFP2", "Jul", 0.51, 0.239, 149], - ["Type I", "Feb", 0.57, 0.272, 162], + ["No mask", "Jul", 10.966992, 12.357222, 809], + ["Type I", "Jul", 2.084648, 1.137819, 149], + ["FFP2", "Jul", 0.654566, .309637, 149], + ["Type I", "Feb", 0.612366, 0.272, 162], ], ) def test_small_shared_office_Geneva(mask_type, month, expected_pi, From 64ebb663e70d0ded43d9c80f747a210609aa8360 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 25 Apr 2022 15:14:46 +0200 Subject: [PATCH 16/38] Handled mypy errors and increased number of decimals on test exposure model --- cara/models.py | 4 ++-- cara/tests/models/test_concentration_model.py | 2 +- cara/tests/models/test_exposure_model.py | 14 +++++++------- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/cara/models.py b/cara/models.py index 1adb810e..65f820e9 100644 --- a/cara/models.py +++ b/cara/models.py @@ -461,8 +461,8 @@ class SARSCoV2(Virus): """ # Updated to use the formula from Dabish et al. https://doi.org/10.1080/02786826.2020.1829536 # with a minimum at hl = 1.1 - inside_temp = inside_temp.value(time) - return max(1.1, (0.693/(0.16030 + 0.04018*(((inside_temp-273.15)-20.615)/10.585)+0.02176*((humidity-45.235)/28.665)+0.1))) + temperature: _VectorisedFloat = inside_temp.value(time) + return np.maximum(1.1, (0.693/(0.16030 + 0.04018*(((temperature-273.15)-20.615)/10.585)+0.02176*((humidity-45.235)/28.665)+0.1))) Virus.types = { diff --git a/cara/tests/models/test_concentration_model.py b/cara/tests/models/test_concentration_model.py index 3b1d9292..ea11a005 100644 --- a/cara/tests/models/test_concentration_model.py +++ b/cara/tests/models/test_concentration_model.py @@ -59,7 +59,7 @@ def test_concentration_model_vectorisation(override_params): def simple_conc_model(): interesting_times = models.SpecificInterval(([0.5, 1.], [1.1, 2], [2., 3.]), ) return models.ConcentrationModel( - models.Room(75), + models.Room(75, models.PiecewiseConstant((0., 24.), (293,))), models.AirChange(interesting_times, 100), models.InfectedPopulation( number=1, diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index abc476d0..428de12d 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -152,7 +152,7 @@ def conc_model(): ) always = models.SpecificInterval(((0., 24.), )) return models.ConcentrationModel( - models.Room(25), + models.Room(25, models.PiecewiseConstant((0., 24.), (293,))), models.AirChange(always, 5), models.EmittingPopulation( number=1, @@ -179,12 +179,12 @@ def sr_model(): @pytest.mark.parametrize( ["exposed_time_interval", "expected_deposited_exposure"], [ - [(0., 1.), 48.19316], - [(1., 1.01), 0.566368], - [(1.01, 1.02), 0.551401], - [(12., 12.01), 0.016278], - [(12., 24.), 691.21381], - [(0., 24.), 750.258043], + [(0., 1.), 48.193159880096644], + [(1., 1.01), 0.5663683904832492], + [(1.01, 1.02), 0.5514013220457682], + [(12., 12.01), 0.016277647772557004], + [(12., 24.), 691.2138099188268], + [(0., 24.), 750.2580429542696], ] ) def test_exposure_model_integral_accuracy(exposed_time_interval, From 9b137200625ce79eb4534a5a799283bb3986c667 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 25 Apr 2022 15:37:55 +0200 Subject: [PATCH 17/38] Fixed bug on frequency and duration on natural ventilation --- cara/apps/calculator/model_generator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index c4593b49..196487bc 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -324,7 +324,7 @@ class FormData: # Initializes a ventilation instance as a window if 'natural_ventilation' is selected, or as a HEPA-filter otherwise if self.ventilation_type == 'natural_ventilation': if self.window_opening_regime == 'windows_open_periodically': - window_interval = models.PeriodicInterval(self.windows_frequency, self.windows_duration, min(self.infected_start, self.exposed_start)) + window_interval = models.PeriodicInterval(self.windows_frequency, self.windows_duration, min(self.infected_start, self.exposed_start)/60) else: window_interval = always_on From 0a18ce7e96182d71097ab7eb06297017430ce8d3 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 25 Apr 2022 15:52:34 +0200 Subject: [PATCH 18/38] fixed test values --- cara/tests/apps/calculator/test_model_generator.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/cara/tests/apps/calculator/test_model_generator.py b/cara/tests/apps/calculator/test_model_generator.py index 5cc9e409..3bfb08f5 100644 --- a/cara/tests/apps/calculator/test_model_generator.py +++ b/cara/tests/apps/calculator/test_model_generator.py @@ -59,7 +59,7 @@ def test_ventilation_slidingwindow(baseline_form: model_generator.FormData): assert isinstance(baseline_window, models.SlidingWindow) window = models.SlidingWindow( - active=models.PeriodicInterval(period=120, duration=10, start=minutes_since_midnight(9 * 60)), + active=models.PeriodicInterval(period=120, duration=10, start=9), inside_temp=models.PiecewiseConstant((0, 24), (293,)), outside_temp=baseline_window.outside_temp, window_height=1.6, opening_length=0.6, @@ -91,7 +91,7 @@ def test_ventilation_hingedwindow(baseline_form: model_generator.FormData): assert isinstance(baseline_window, models.HingedWindow) window = models.HingedWindow( - active=models.PeriodicInterval(period=120, duration=10, start=minutes_since_midnight(9 * 60)), + active=models.PeriodicInterval(period=120, duration=10, start=9), inside_temp=models.PiecewiseConstant((0, 24), (293,)), outside_temp=baseline_window.outside_temp, window_height=1.6, window_width=1., opening_length=0.6, @@ -152,7 +152,7 @@ def test_ventilation_window_hepa(baseline_form: model_generator.FormData): # Now build the equivalent ventilation instance directly, and compare. window = models.SlidingWindow( - active=models.PeriodicInterval(period=120, duration=10, start=minutes_since_midnight(9 * 60)), + active=models.PeriodicInterval(period=120, duration=10, start=9), inside_temp=models.PiecewiseConstant((0, 24), (293,)), outside_temp=baseline_window.outside_temp, window_height=1.6, opening_length=0.6, From 71078671f99dee879bfbcf36ff033466feeb26f5 Mon Sep 17 00:00:00 2001 From: jdevine Date: Mon, 25 Apr 2022 16:20:29 +0200 Subject: [PATCH 19/38] Updated formula to reflect published correction --- cara/models.py | 7 +++++-- cara/tests/test_full_algorithm.py | 7 ++++--- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/cara/models.py b/cara/models.py index 65f820e9..f2ac2500 100644 --- a/cara/models.py +++ b/cara/models.py @@ -459,10 +459,13 @@ class SARSCoV2(Virus): piecewise constant model (for more details see A. Henriques et al, CERN-OPEN-2021-004, DOI: 10.17181/CERN.1GDQ.5Y75) """ - # Updated to use the formula from Dabish et al. https://doi.org/10.1080/02786826.2020.1829536 + # Updated to use the formula from Dabish et al. with correction https://doi.org/10.1080/02786826.2020.1829536 # with a minimum at hl = 1.1 temperature: _VectorisedFloat = inside_temp.value(time) - return np.maximum(1.1, (0.693/(0.16030 + 0.04018*(((temperature-273.15)-20.615)/10.585)+0.02176*((humidity-45.235)/28.665)+0.1))) + return np.maximum(1.1, (0.693/((0.16030 + 0.04018*(((temperature-273.15)-20.615)/10.585) + +0.02176*((humidity-45.235)/28.665) + -0.14369 + -0.2636*((temperature-273.15)-20.615)/10.585)))) Virus.types = { diff --git a/cara/tests/test_full_algorithm.py b/cara/tests/test_full_algorithm.py index 71f10f05..8a07321c 100644 --- a/cara/tests/test_full_algorithm.py +++ b/cara/tests/test_full_algorithm.py @@ -89,9 +89,10 @@ class SimpleConcentrationModel: """ return (self.lambda_ventilation - + ln2/(max(1.1, (0.693 / (0.16030 + 0.04018 * (((22) - 20.615) / 10.585) + 0.02176 * ( - (self.humidity - 45.235) / 28.665) + 0.1))))) - #6.43 if self.humidity<=0.4 else 1.1) ) + + ln2/(max(1.1, (0.693 / ((0.16030 + 0.04018 * (((22) - 20.615) / 10.585) + + 0.02176 * ((self.humidity - 45.235) / 28.665) + - 0.14369 + - 0.2636((22-20.615)/10.585))))))) @method_cache def deposition_removal_coefficient(self) -> float: From ac49e102d4bb3e1ce02673275b8a7a4b11b29c6d Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 25 Apr 2022 15:44:27 +0200 Subject: [PATCH 20/38] Removing commented decorator in test_full_algorithm.py; inverting order Speaking/Breathing in short-range tests in test_full_algorithm.py --- cara/tests/test_full_algorithm.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/cara/tests/test_full_algorithm.py b/cara/tests/test_full_algorithm.py index 6726f3a4..372ca8b1 100644 --- a/cara/tests/test_full_algorithm.py +++ b/cara/tests/test_full_algorithm.py @@ -19,7 +19,7 @@ from cara.monte_carlo.data import (expiration_distributions, # TODO: seed better the random number generators np.random.seed(2000) SAMPLE_SIZE = 1_000_000 -TOLERANCE = 0.02 +TOLERANCE = 0.04 sqrt2pi = np.sqrt(2.*np.pi) sqrt2 = np.sqrt(2.) @@ -244,7 +244,6 @@ class SimpleShortRangeModel: return dilution -# @method_cache def jet_concentration(self,conc_model: SimpleConcentrationModel) -> _VectorisedFloat: """ virion concentration at the origin of the jet (close to @@ -500,14 +499,14 @@ def c_model_distr() -> mc.ConcentrationModel: def sr_models() -> typing.Tuple[mc.ShortRangeModel, ...]: return ( mc.ShortRangeModel( - expiration = short_range_expiration_distributions['Breathing'], + expiration = short_range_expiration_distributions['Speaking'], 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'], + expiration = short_range_expiration_distributions['Breathing'], + activity = models.Activity.types['Heavy exercise'], presence = interaction_intervals[1], distance = 0.854, ).build_model(SAMPLE_SIZE), @@ -533,14 +532,14 @@ def simple_sr_models() -> typing.Tuple[SimpleShortRangeModel, ...]: interaction_interval = interaction_intervals[0], distance = 0.854, breathing_rate = models.Activity.types['Seated'].exhalation_rate, - BLO_factors = expiration_BLO_factors['Breathing'], + BLO_factors = expiration_BLO_factors['Speaking'], ), SimpleShortRangeModel( interaction_interval = interaction_intervals[1], distance = 0.854, - breathing_rate = models.Activity.types['Seated'].exhalation_rate, - BLO_factors = expiration_BLO_factors['Speaking'], - ) + breathing_rate = models.Activity.types['Heavy exercise'].exhalation_rate, + BLO_factors = expiration_BLO_factors['Breathing'], + ), ) From c665f31df98ff4089eb0b0c03974cf2bc9d87b7d Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 25 Apr 2022 17:24:03 +0200 Subject: [PATCH 21/38] adjusted tests for the formula from Dabish et al --- cara/models.py | 2 +- cara/tests/models/test_exposure_model.py | 12 ++++++------ cara/tests/test_full_algorithm.py | 11 ++++------- cara/tests/test_known_quantities.py | 12 ++++++------ cara/tests/test_monte_carlo_full_models.py | 20 ++++++++++---------- 5 files changed, 27 insertions(+), 30 deletions(-) diff --git a/cara/models.py b/cara/models.py index f2ac2500..0a94a2dd 100644 --- a/cara/models.py +++ b/cara/models.py @@ -465,7 +465,7 @@ class SARSCoV2(Virus): return np.maximum(1.1, (0.693/((0.16030 + 0.04018*(((temperature-273.15)-20.615)/10.585) +0.02176*((humidity-45.235)/28.665) -0.14369 - -0.2636*((temperature-273.15)-20.615)/10.585)))) + -0.02636*((temperature-273.15)-20.615)/10.585)))) Virus.types = { diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 428de12d..5a6a5449 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -179,12 +179,12 @@ def sr_model(): @pytest.mark.parametrize( ["exposed_time_interval", "expected_deposited_exposure"], [ - [(0., 1.), 48.193159880096644], - [(1., 1.01), 0.5663683904832492], - [(1.01, 1.02), 0.5514013220457682], - [(12., 12.01), 0.016277647772557004], - [(12., 24.), 691.2138099188268], - [(0., 24.), 750.2580429542696], + [(0., 1.), 45.6008710], + [(1., 1.01), 0.5280401], + [(1.01, 1.02), 0.51314096385], + [(12., 12.01), 0.016255813185], + [(12., 24.), 645.63619275], + [(0., 24.), 700.7322474], ] ) def test_exposure_model_integral_accuracy(exposed_time_interval, diff --git a/cara/tests/test_full_algorithm.py b/cara/tests/test_full_algorithm.py index 8a07321c..3f701e5e 100644 --- a/cara/tests/test_full_algorithm.py +++ b/cara/tests/test_full_algorithm.py @@ -47,9 +47,6 @@ class SimpleConcentrationModel: #: room volume (m^3) room_volume: _VectorisedFloat - #: The temperature inside the room (Kelvin). - #inside_temp: float = 293.15 - #: ventilation rate (air changes per hour) - including HEPA lambda_ventilation: _VectorisedFloat @@ -89,10 +86,10 @@ class SimpleConcentrationModel: """ return (self.lambda_ventilation - + ln2/(max(1.1, (0.693 / ((0.16030 + 0.04018 * (((22) - 20.615) / 10.585) - + 0.02176 * ((self.humidity - 45.235) / 28.665) + + ln2/(np.maximum(1.1, (0.693 / ((0.16030 + 0.04018 * (((21) - 20.615) / 10.585) + + 0.02176*((self.humidity - 45.235) / 28.665) - 0.14369 - - 0.2636((22-20.615)/10.585))))))) + - 0.02636*((21-20.615)/10.585))))))) @method_cache def deposition_removal_coefficient(self) -> float: @@ -460,7 +457,7 @@ interaction_intervals = (models.SpecificInterval(present_times=((10.5, 11.0),)), @pytest.fixture def c_model() -> mc.ConcentrationModel: return mc.ConcentrationModel( - room=models.Room(volume=50, humidity=0.3), + room=models.Room(volume=50, inside_temp=models.PiecewiseConstant((0., 24.), (293,)), humidity=0.3), ventilation=models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=1.), infected=mc.InfectedPopulation( number=1, diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 4eb94d0b..ed775312 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -43,7 +43,7 @@ def test_concentrations(baseline_concentration_model): concentrations = [baseline_concentration_model.concentration(float(t)) for t in ts] npt.assert_allclose( concentrations, - [0.000000e+00, 2.108144e+01, 1.004741e-12, 2.108144e+01, 4.788592e-26], + [0.000000e+00, 20.805628, 6.602814e-13, 20.805628, 2.09545e-26], rtol=1e-6 ) @@ -94,7 +94,7 @@ def test_r0(baseline_exposure_model): # expected r0 was computed with a trapezoidal integration, using # a mesh of 100'000 pts per exposed presence interval. r0 = baseline_exposure_model.reproduction_number() - npt.assert_allclose(r0, 781.293793) + npt.assert_allclose(r0, 776.941990) def test_periodic_window(baseline_periodic_window, baseline_room): @@ -381,8 +381,8 @@ def build_exposure_model(concentration_model, short_range_model): @pytest.mark.parametrize( "month, expected_deposited_exposure", [ - ['Jan', 392.994454], - ['Jun', 2127.82386], + ['Jan', 377.440565819], + ['Jun', 1721.03336729], ], ) def test_exposure_hourly_dep(month,expected_deposited_exposure, baseline_sr_model): @@ -402,8 +402,8 @@ def test_exposure_hourly_dep(month,expected_deposited_exposure, baseline_sr_mode @pytest.mark.parametrize( "month, expected_deposited_exposure", [ - ['Jan', 394.000261], - ['Jun', 2239.777906], + ['Jan', 383.339206111], + ['Jun', 1799.17597184], ], ) def test_exposure_hourly_dep_refined(month,expected_deposited_exposure, baseline_sr_model): diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index 58f0a57d..48682758 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -310,13 +310,13 @@ def waiting_room_mc(): @pytest.mark.parametrize( "mc_model, expected_pi, expected_new_cases, expected_dose, expected_ER", [ - ["shared_office_mc", 6.524158, 0.195725, 3.413983, 809], - ["classroom_mc", 9.5, 1.85, 9.478, 5624], - ["ski_cabin_mc", 16.0, 0.5, 17.315, 7966], - ["skagit_chorale_mc",69.84901, 40.0, 121.265911, 190422], - ["bus_ride_mc", 13.311721, 8.918853, 8.6662, 5419], + ["shared_office_mc", 6.03, 0.18, 3.198, 809], + ["classroom_mc", 8.63, 1.64, 8.082, 5624], + ["ski_cabin_mc", 16.0, 0.47, 17.315, 7966], + ["skagit_chorale_mc",65.7, 40.0, 102.213, 190422], + ["bus_ride_mc", 12.0, 8.0, 7.65, 5419], ["gym_mc", 0.45, 0.13, 0.208, 1145], - ["waiting_room_mc", 1.854373, 0.259612, 0.99173, 737], + ["waiting_room_mc", 1.59, 0.22, 0.821, 737], ] ) def test_report_models(mc_model, expected_pi, expected_new_cases, @@ -337,10 +337,10 @@ def test_report_models(mc_model, expected_pi, expected_new_cases, @pytest.mark.parametrize( "mask_type, month, expected_pi, expected_dose, expected_ER", [ - ["No mask", "Jul", 10.966992, 12.357222, 809], - ["Type I", "Jul", 2.084648, 1.137819, 149], - ["FFP2", "Jul", 0.654566, .309637, 149], - ["Type I", "Feb", 0.612366, 0.272, 162], + ["No mask", "Jul", 9.52, 9.920, 809], + ["Type I", "Jul", 1.7, 0.913, 149], + ["FFP2", "Jul", 0.51, 0.239, 149], + ["Type I", "Feb", 0.57, 0.272, 162], ], ) def test_small_shared_office_Geneva(mask_type, month, expected_pi, From 2d1f7fa823f79f6a58f71cedd33b81ab814e9868 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 25 Apr 2022 17:25:59 +0200 Subject: [PATCH 22/38] Fixing bug pointed in issue #262 --- cara/models.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/cara/models.py b/cara/models.py index cbb547c9..2c65c770 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1356,7 +1356,7 @@ class ExposureModel: # to perform properly the Monte-Carlo integration over # particle diameters (doing things in another order would # lead to wrong results for the probability of infection). - deposited_exposure += (np.array(short_range_jet_exposure + this_deposited_exposure = (np.array(short_range_jet_exposure * fdep).mean()/dilution - np.array(short_range_lr_exposure * fdep).mean() * self.concentration_model.infected.activity.exhalation_rate @@ -1364,14 +1364,15 @@ class ExposureModel: 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_jet_exposure + this_deposited_exposure = (short_range_jet_exposure * fdep/dilution - short_range_lr_exposure * fdep * self.concentration_model.infected.activity.exhalation_rate /dilution) # multiply by the (diameter-independent) inhalation rate - deposited_exposure *= interaction.activity.inhalation_rate + deposited_exposure += (this_deposited_exposure * + interaction.activity.inhalation_rate) # then we multiply by diameter-independent quantities: viral load # and fraction of infected virions From faaf1ed086068ee9f845fa289977bea0f9c08e11 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Mon, 25 Apr 2022 17:32:30 +0200 Subject: [PATCH 23/38] In models, ExposureModel: multiplication by dilution factor done in a better place, in deposited_exposure_between_bounds method --- cara/models.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/cara/models.py b/cara/models.py index 2c65c770..9d7898fd 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1357,22 +1357,20 @@ class ExposureModel: # particle diameters (doing things in another order would # lead to wrong results for the probability of infection). this_deposited_exposure = (np.array(short_range_jet_exposure - * fdep).mean()/dilution + * fdep).mean() - np.array(short_range_lr_exposure * fdep).mean() - * self.concentration_model.infected.activity.exhalation_rate - /dilution) + * self.concentration_model.infected.activity.exhalation_rate) else: # in the case of a single diameter or no diameter defined, # one should not take any mean at this stage. - this_deposited_exposure = (short_range_jet_exposure - * fdep/dilution + this_deposited_exposure = (short_range_jet_exposure * fdep - short_range_lr_exposure * fdep - * self.concentration_model.infected.activity.exhalation_rate - /dilution) + * self.concentration_model.infected.activity.exhalation_rate) # multiply by the (diameter-independent) inhalation rate deposited_exposure += (this_deposited_exposure * - interaction.activity.inhalation_rate) + interaction.activity.inhalation_rate + /dilution) # then we multiply by diameter-independent quantities: viral load # and fraction of infected virions From f0bdda7c7d991e35ad705c424b1ae68cca7728c1 Mon Sep 17 00:00:00 2001 From: Nicolas Mounet Date: Tue, 26 Apr 2022 11:19:33 +0200 Subject: [PATCH 24/38] Adding mask inhale efficiency factor to short-range exposure; adding corresponding test; some cosmetics in test_short_range_model.py --- cara/models.py | 2 +- cara/tests/models/test_short_range_model.py | 58 ++++++++++++++++++--- 2 files changed, 51 insertions(+), 9 deletions(-) diff --git a/cara/models.py b/cara/models.py index 9d7898fd..c6df3464 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1377,7 +1377,7 @@ class ExposureModel: f_inf = self.concentration_model.infected.fraction_of_infectious_virus() deposited_exposure *= (f_inf * self.concentration_model.virus.viral_load_in_sputum - ) + * (1 - self.exposed.mask.inhale_efficiency())) # long-range concentration deposited_exposure += self.long_range_deposited_exposure_between_bounds(time1, time2) diff --git a/cara/tests/models/test_short_range_model.py b/cara/tests/models/test_short_range_model.py index b1dee050..ceef0bb2 100644 --- a/cara/tests/models/test_short_range_model.py +++ b/cara/tests/models/test_short_range_model.py @@ -6,10 +6,13 @@ 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 +from cara.monte_carlo.data import short_range_expiration_distributions,\ + expiration_distributions, short_range_distances, activity_distributions # TODO: seed better the random number generators np.random.seed(2000) +SAMPLE_SIZE = 250_000 + @pytest.fixture def concentration_model() -> mc_models.ConcentrationModel: @@ -41,8 +44,8 @@ def short_range_model(): 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) + concentration_model = concentration_model.build_model(SAMPLE_SIZE) + model = short_range_model.build_model(SAMPLE_SIZE) 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_jet_exposure_between_bounds(concentration_model, 10.75, 10.85), np.ndarray) @@ -60,10 +63,10 @@ def test_short_range_model_ndarray(concentration_model, short_range_model): ] ) def test_dilution_factor(activity, expected_dilution): - model = models.ShortRangeModel(expiration="Breathing", + model = mc_models.ShortRangeModel(expiration=short_range_expiration_distributions['Breathing'], activity=models.Activity.types[activity], presence=models.SpecificInterval(present_times=((10.5, 11.0),)), - distance=0.854) + distance=0.854).build_model(SAMPLE_SIZE) assert isinstance(model.dilution_factor(), np.ndarray) np.testing.assert_almost_equal( model.dilution_factor(), expected_dilution, decimal=10 @@ -105,11 +108,50 @@ def test_extract_between_bounds(short_range_model, time1, time2, [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) +def test_short_range_concentration(time, expected_short_range_concentration, + concentration_model, short_range_model): + concentration_model = concentration_model.build_model(SAMPLE_SIZE) + model = short_range_model.build_model(SAMPLE_SIZE) np.testing.assert_allclose( np.array(model.short_range_concentration(concentration_model, time)).mean(), expected_short_range_concentration, rtol=0.02 ) + +def test_short_range_exposure_with_ndarray_mask(): + c_model = mc_models.ConcentrationModel( + room=models.Room(volume=50, humidity=0.3), + ventilation=models.AirChange(active=models.PeriodicInterval(period=120, duration=120), + air_exch=10_000_000,), + infected=mc_models.InfectedPopulation( + number=1, + presence=models.SpecificInterval(present_times=((8.5, 12.5), (13.5, 17.5))), + 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, + ) + sr_model = mc_models.ShortRangeModel(expiration=short_range_expiration_distributions['Shouting'], + activity=models.Activity.types['Heavy exercise'], + presence=models.SpecificInterval(present_times=((10.5, 11.0),)), + distance=0.854) + e_model = mc_models.ExposureModel( + concentration_model = c_model, + short_range = (sr_model,), + exposed = mc_models.Population( + number=1, + presence=models.SpecificInterval(present_times=((8.5, 12.5), (13.5, 17.5))), + mask=models.Mask(η_inhale=np.array([0., 0.3, 0.5])), + activity=models.Activity.types['Light activity'], + host_immunity=0., + ), + ).build_model(SAMPLE_SIZE) + assert isinstance(e_model.deposited_exposure(), np.ndarray) + assert len(e_model.deposited_exposure()) == 3 + np.testing.assert_allclose(e_model.deposited_exposure(), + e_model.deposited_exposure()[0]*np.array([1., 0.7, 0.5]), + rtol=1e-8) + From 71f36e0aa41d5acab314ab5d2a2f03fb62e5845b Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 27 Apr 2022 14:09:24 +0200 Subject: [PATCH 25/38] Changed inside_temp to be passed as vectorised_float in the decay_constant method call --- cara/models.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/cara/models.py b/cara/models.py index 0a94a2dd..cb389be0 100644 --- a/cara/models.py +++ b/cara/models.py @@ -294,7 +294,6 @@ class WindowOpening(Ventilation): def transition_times(self, room: Room) -> typing.Set[float]: transitions = super().transition_times(room) - print(room.inside_temp.transition_times) transitions.update(room.inside_temp.transition_times) transitions.update(self.outside_temp.transition_times) return transitions @@ -440,20 +439,20 @@ class Virus: #: Pre-populated examples of Viruses. types: typing.ClassVar[typing.Dict[str, "Virus"]] - def halflife(self, humidity: _VectorisedFloat, inside_temp: PiecewiseConstant, time: float) -> _VectorisedFloat: + def halflife(self, humidity: _VectorisedFloat, inside_temp: _VectorisedFloat) -> _VectorisedFloat: # Biological decay (inactivation of the virus in air) - virus # dependent and function of humidity raise NotImplementedError - def decay_constant(self, humidity: _VectorisedFloat, inside_temp: PiecewiseConstant, time: float) -> _VectorisedFloat: + def decay_constant(self, humidity: _VectorisedFloat, inside_temp: _VectorisedFloat) -> _VectorisedFloat: # Viral inactivation per hour (h^-1) (function of humidity) - return np.log(2) / self.halflife(humidity, inside_temp, time) + return np.log(2) / self.halflife(humidity, inside_temp) @dataclass(frozen=True) class SARSCoV2(Virus): - def halflife(self, humidity: _VectorisedFloat, inside_temp: PiecewiseConstant, time: float) -> _VectorisedFloat: + def halflife(self, humidity: _VectorisedFloat, inside_temp: _VectorisedFloat) -> _VectorisedFloat: """ Half-life changes with humidity level. Here is implemented a simple piecewise constant model (for more details see A. Henriques et al, @@ -461,11 +460,10 @@ class SARSCoV2(Virus): """ # Updated to use the formula from Dabish et al. with correction https://doi.org/10.1080/02786826.2020.1829536 # with a minimum at hl = 1.1 - temperature: _VectorisedFloat = inside_temp.value(time) - return np.maximum(1.1, (0.693/((0.16030 + 0.04018*(((temperature-273.15)-20.615)/10.585) + return np.maximum(1.1, (0.693/((0.16030 + 0.04018*(((inside_temp-273.15)-20.615)/10.585) +0.02176*((humidity-45.235)/28.665) -0.14369 - -0.02636*((temperature-273.15)-20.615)/10.585)))) + -0.02636*((inside_temp-273.15)-20.615)/10.585)))) Virus.types = { @@ -924,7 +922,7 @@ class ConcentrationModel: k = (vg * 3600) / h #todo: Inside_temp needs to be exposed/added to the room; return ( - k + self.virus.decay_constant(self.room.humidity, self.room.inside_temp, time) + k + self.virus.decay_constant(self.room.humidity, self.room.inside_temp.value(time)) + self.ventilation.air_exchange(self.room, time) ) From 5dcb698622666c293692e7c986b9e4c74c6d6d4f Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 27 Apr 2022 15:23:46 +0200 Subject: [PATCH 26/38] Updated humidity value to be in percentage --- cara/models.py | 6 +++--- cara/tests/models/test_exposure_model.py | 12 +++++------ cara/tests/models/test_virus.py | 24 ++++++++++++++++++++++ cara/tests/test_full_algorithm.py | 2 +- cara/tests/test_known_quantities.py | 12 +++++------ cara/tests/test_monte_carlo_full_models.py | 20 +++++++++--------- 6 files changed, 50 insertions(+), 26 deletions(-) create mode 100644 cara/tests/models/test_virus.py diff --git a/cara/models.py b/cara/models.py index cb389be0..333f84bf 100644 --- a/cara/models.py +++ b/cara/models.py @@ -445,7 +445,7 @@ class Virus: raise NotImplementedError def decay_constant(self, humidity: _VectorisedFloat, inside_temp: _VectorisedFloat) -> _VectorisedFloat: - # Viral inactivation per hour (h^-1) (function of humidity) + # Viral inactivation per hour (h^-1) (function of humidity and inside temperature) return np.log(2) / self.halflife(humidity, inside_temp) @@ -459,9 +459,9 @@ class SARSCoV2(Virus): CERN-OPEN-2021-004, DOI: 10.17181/CERN.1GDQ.5Y75) """ # Updated to use the formula from Dabish et al. with correction https://doi.org/10.1080/02786826.2020.1829536 - # with a minimum at hl = 1.1 + # with a minimum at hl = 1.1. Note that humidity is in percentage and inside_temp in °C. return np.maximum(1.1, (0.693/((0.16030 + 0.04018*(((inside_temp-273.15)-20.615)/10.585) - +0.02176*((humidity-45.235)/28.665) + +0.02176*(((humidity*100)-45.235)/28.665) -0.14369 -0.02636*((inside_temp-273.15)-20.615)/10.585)))) diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 5a6a5449..b1959ca8 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -179,12 +179,12 @@ def sr_model(): @pytest.mark.parametrize( ["exposed_time_interval", "expected_deposited_exposure"], [ - [(0., 1.), 45.6008710], - [(1., 1.01), 0.5280401], - [(1.01, 1.02), 0.51314096385], - [(12., 12.01), 0.016255813185], - [(12., 24.), 645.63619275], - [(0., 24.), 700.7322474], + [(0., 1.), 49.60355486823376], + [(1., 1.01), 0.5876509666617122], + [(1.01, 1.02), 0.5726560233646302], + [(12., 12.01), 0.016288631412456556], + [(12., 24.), 716.6229828782525], + [(0., 24.), 777.8717785392312], ] ) def test_exposure_model_integral_accuracy(exposed_time_interval, diff --git a/cara/tests/models/test_virus.py b/cara/tests/models/test_virus.py new file mode 100644 index 00000000..0f875fa8 --- /dev/null +++ b/cara/tests/models/test_virus.py @@ -0,0 +1,24 @@ +import numpy as np +import numpy.testing as npt +import pytest + +from cara import models + +@pytest.mark.parametrize( + "inside_temp, humidity, expected_halflife, expected_decay_constant", + [ + [293.15, 0.5, 35.67710693238622, 0.01942834607844098], + [272.15, 0.7, 96.40459058258793, 0.007189981061805762], + [300.15, 1., 10.418034697539541, 0.0665333914393324], + [300.15, 0., 1.1, 0.6301338005090411], + [np.array([272.15, 300.15]), np.array([0.7, 0.]), + np.array([96.40459058258793, 1.1]), np.array([0.007189981061805762, 0.6301338005090411])], + [np.array([293.15, 300.15]), np.array([0.5, 1.]), + np.array([35.67710693238622, 10.418034697539541]), np.array([0.01942834607844098, 0.0665333914393324])] + ], +) +def test_decay_constant(inside_temp, humidity, expected_halflife, expected_decay_constant): + npt.assert_equal(models.Virus.types['SARS_CoV_2'].halflife(humidity, inside_temp), + expected_halflife) + npt.assert_equal(models.Virus.types['SARS_CoV_2'].decay_constant(humidity, inside_temp), + expected_decay_constant) \ No newline at end of file diff --git a/cara/tests/test_full_algorithm.py b/cara/tests/test_full_algorithm.py index 3f701e5e..1fac1c8b 100644 --- a/cara/tests/test_full_algorithm.py +++ b/cara/tests/test_full_algorithm.py @@ -87,7 +87,7 @@ class SimpleConcentrationModel: return (self.lambda_ventilation + ln2/(np.maximum(1.1, (0.693 / ((0.16030 + 0.04018 * (((21) - 20.615) / 10.585) - + 0.02176*((self.humidity - 45.235) / 28.665) + + 0.02176*(((self.humidity * 100) - 45.235) / 28.665) - 0.14369 - 0.02636*((21-20.615)/10.585))))))) diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index ed775312..98a90acb 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -43,7 +43,7 @@ def test_concentrations(baseline_concentration_model): 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], + [0.000000e+00, 2.122276e+01, 1.240684e-12, 2.122276e+01, 7.253047e-26], rtol=1e-6 ) @@ -94,7 +94,7 @@ def test_r0(baseline_exposure_model): # expected r0 was computed with a trapezoidal integration, using # a mesh of 100'000 pts per exposed presence interval. r0 = baseline_exposure_model.reproduction_number() - npt.assert_allclose(r0, 776.941990) + npt.assert_allclose(r0, 783.490035) def test_periodic_window(baseline_periodic_window, baseline_room): @@ -381,8 +381,8 @@ def build_exposure_model(concentration_model, short_range_model): @pytest.mark.parametrize( "month, expected_deposited_exposure", [ - ['Jan', 377.440565819], - ['Jun', 1721.03336729], + ['Jan', 401.300989], + ['Jun', 2420.383151], ], ) def test_exposure_hourly_dep(month,expected_deposited_exposure, baseline_sr_model): @@ -402,8 +402,8 @@ def test_exposure_hourly_dep(month,expected_deposited_exposure, baseline_sr_mode @pytest.mark.parametrize( "month, expected_deposited_exposure", [ - ['Jan', 383.339206111], - ['Jun', 1799.17597184], + ['Jan', 402.348745], + ['Jun', 2558.632473], ], ) def test_exposure_hourly_dep_refined(month,expected_deposited_exposure, baseline_sr_model): diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index 48682758..d60dbcaf 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -310,13 +310,13 @@ def waiting_room_mc(): @pytest.mark.parametrize( "mc_model, expected_pi, expected_new_cases, expected_dose, expected_ER", [ - ["shared_office_mc", 6.03, 0.18, 3.198, 809], - ["classroom_mc", 8.63, 1.64, 8.082, 5624], + ["shared_office_mc", 6.75, 0.20, 3.594, 809], + ["classroom_mc", 9.58, 1.82, 9.664, 5624], ["ski_cabin_mc", 16.0, 0.47, 17.315, 7966], - ["skagit_chorale_mc",65.7, 40.0, 102.213, 190422], - ["bus_ride_mc", 12.0, 8.0, 7.65, 5419], - ["gym_mc", 0.45, 0.13, 0.208, 1145], - ["waiting_room_mc", 1.59, 0.22, 0.821, 737], + ["skagit_chorale_mc",72.20, 43.32, 134.234, 190422], + ["bus_ride_mc", 14.08, 9.43, 9.26, 5419], + ["gym_mc", 0.49, 0.14, 0.226, 1145], + ["waiting_room_mc", 2.02, 0.28, 1.104, 737], ] ) def test_report_models(mc_model, expected_pi, expected_new_cases, @@ -337,10 +337,10 @@ def test_report_models(mc_model, expected_pi, expected_new_cases, @pytest.mark.parametrize( "mask_type, month, expected_pi, expected_dose, expected_ER", [ - ["No mask", "Jul", 9.52, 9.920, 809], - ["Type I", "Jul", 1.7, 0.913, 149], - ["FFP2", "Jul", 0.51, 0.239, 149], - ["Type I", "Feb", 0.57, 0.272, 162], + ["No mask", "Jul", 11.81, 14.137, 809], + ["Type I", "Jul", 2.33, 1.305, 149], + ["FFP2", "Jul", 0.73, 0.351, 149], + ["Type I", "Feb", 0.62, 0.291, 162], ], ) def test_small_shared_office_Geneva(mask_type, month, expected_pi, From 236eaee456b9aa57b0e5261d2c0aec75dfdd7b55 Mon Sep 17 00:00:00 2001 From: Nicola Tarocco Date: Thu, 28 Apr 2022 15:44:02 +0200 Subject: [PATCH 27/38] tests: increase request timeout * increases HTTP request timeout in tests * closes #257 --- cara/tests/apps/calculator/test_webapp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cara/tests/apps/calculator/test_webapp.py b/cara/tests/apps/calculator/test_webapp.py index 65f3b216..e6b4eedf 100644 --- a/cara/tests/apps/calculator/test_webapp.py +++ b/cara/tests/apps/calculator/test_webapp.py @@ -6,7 +6,7 @@ import tornado.testing import cara.apps.calculator from cara.apps.calculator.report_generator import generate_permalink -_TIMEOUT = 20. +_TIMEOUT = 30. @pytest.fixture def app(): From a7bcc3ec2cb33d183bd427aa9053b94a64d79814 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 2 May 2022 11:18:49 +0200 Subject: [PATCH 28/38] Removed beta distribution and added custom distribution for the interpersonal distances distribution --- cara/monte_carlo/data.py | 8 ++++++-- cara/monte_carlo/sampleable.py | 19 ------------------- cara/tests/models/test_short_range_model.py | 6 +++--- 3 files changed, 9 insertions(+), 24 deletions(-) diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index a725b8fa..bf9c148c 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -5,7 +5,7 @@ import numpy as np from scipy import special as sp import cara.monte_carlo as mc -from cara.monte_carlo.sampleable import LogNormal,LogCustomKernel,CustomKernel,Uniform, Beta +from cara.monte_carlo.sampleable import LogNormal,LogCustomKernel,CustomKernel,Uniform, Custom sqrt2pi = np.sqrt(2.*np.pi) @@ -203,4 +203,8 @@ short_range_expiration_distributions = { # Derived from Fig 8 a) "stand-stand" in https://www.mdpi.com/1660-4601/17/4/1445/htm -short_range_distances = Beta(alpha=0.588715, beta=1.50214, loc=0.5, scale=1.5) \ No newline at end of file +distances = np.array((0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2)) +frequencies = np.array((0.0598036,0.0946154,0.1299152,0.1064905,0.1099066,0.0998209, 0.0845298,0.0479286,0.0406084,0.039795,0.0205997,0.0152316,0.0118155,0.0118155,0.018485,0.0205997)) +short_range_distances = Custom(bounds=(0.5,2.), + function=lambda x: np.interp(x,distances,frequencies,left=0.,right=0.), + max_function=0.13) \ No newline at end of file diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index d9d4afc7..27907e49 100644 --- a/cara/monte_carlo/sampleable.py +++ b/cara/monte_carlo/sampleable.py @@ -2,7 +2,6 @@ import typing import numpy as np from sklearn.neighbors import KernelDensity # type: ignore -from scipy.stats import beta import cara.models @@ -131,24 +130,6 @@ class LogCustomKernel(SampleableDistribution): return 10 ** kde_model.sample(n_samples=size)[:, 0] -class Beta(SampleableDistribution): - """ - Defines a Beta distribution parameterized by two positive shape parameters, - denoted by alpha (α) and beta (β), that appear as exponents of the random - variable and control the shape of the distribution. - """ - - def __init__(self, alpha: float, beta: float, loc: float, scale: float): - # these are resp. the alpha and beta of the underlying distribution - self.alpha = alpha - self.beta = beta - self.loc = loc - self.scale = scale - - def generate_samples(self, size: int) -> float_array_size_n: - return beta.rvs(a = self.alpha, b = self.beta, loc=self.loc, scale=self.scale, size=size) - - _VectorisedFloatOrSampleable = typing.Union[ SampleableDistribution, cara.models._VectorisedFloat, ] diff --git a/cara/tests/models/test_short_range_model.py b/cara/tests/models/test_short_range_model.py index 481860a8..6f4de3d2 100644 --- a/cara/tests/models/test_short_range_model.py +++ b/cara/tests/models/test_short_range_model.py @@ -72,9 +72,9 @@ def test_dilution_factor(activity, expected_dilution): @pytest.mark.parametrize( "time, expected_short_range_concentration", [ [8.5, 0.], - [10.5, 8.037883241318065], - [10.6, 8.037883241318065], - [11.0, 8.037883241318065], + [10.5, 5.401601371244907], + [10.6, 5.401601371244907], + [11.0, 5.401601371244907], [12.0, 0.], ] ) From f518a893d506058968c1cd722a6ea5e7a04fbc05 Mon Sep 17 00:00:00 2001 From: jdevine Date: Mon, 9 May 2022 14:04:27 +0200 Subject: [PATCH 29/38] Corrected handling of min-1/hr-1 half life Changed half life 'limit' value to 6.43 (as paper model, used when formula gives negative values). --- cara/models.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/cara/models.py b/cara/models.py index 06f45d25..faee2e3f 100644 --- a/cara/models.py +++ b/cara/models.py @@ -459,11 +459,14 @@ class SARSCoV2(Virus): CERN-OPEN-2021-004, DOI: 10.17181/CERN.1GDQ.5Y75) """ # Updated to use the formula from Dabish et al. with correction https://doi.org/10.1080/02786826.2020.1829536 - # with a minimum at hl = 1.1. Note that humidity is in percentage and inside_temp in °C. - return np.maximum(1.1, (0.693/((0.16030 + 0.04018*(((inside_temp-273.15)-20.615)/10.585) + # with a maximum at hl = 6.43. Note that humidity is in percentage and inside_temp in °C. + hl_calc = ((0.693/((0.16030 + 0.04018*(((inside_temp-273.15)-20.615)/10.585) +0.02176*(((humidity*100)-45.235)/28.665) -0.14369 - -0.02636*((inside_temp-273.15)-20.615)/10.585)))) + -0.02636*((inside_temp-273.15)-20.615)/10.585)))/60) + if (hl_calc <= 0): + hl_calc = 6.43 + return hl_calc Virus.types = { From 1e689af19c8a34acb1aedd5d56165ee8ece62a96 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 9 May 2022 16:29:52 +0200 Subject: [PATCH 30/38] updated half-life formula and respective tests --- cara/models.py | 11 ++++++----- cara/tests/models/test_exposure_model.py | 12 ++++++------ cara/tests/models/test_virus.py | 20 ++++++++++---------- cara/tests/test_full_algorithm.py | 9 +++++---- cara/tests/test_known_quantities.py | 12 ++++++------ cara/tests/test_monte_carlo_full_models.py | 20 ++++++++++---------- 6 files changed, 43 insertions(+), 41 deletions(-) diff --git a/cara/models.py b/cara/models.py index faee2e3f..34c0d666 100644 --- a/cara/models.py +++ b/cara/models.py @@ -459,14 +459,15 @@ class SARSCoV2(Virus): CERN-OPEN-2021-004, DOI: 10.17181/CERN.1GDQ.5Y75) """ # Updated to use the formula from Dabish et al. with correction https://doi.org/10.1080/02786826.2020.1829536 - # with a maximum at hl = 6.43. Note that humidity is in percentage and inside_temp in °C. - hl_calc = ((0.693/((0.16030 + 0.04018*(((inside_temp-273.15)-20.615)/10.585) + # with a maximum at hl = 6.43 (compensate for the negative decay values in the paper). + # Note that humidity is in percentage and inside_temp in °C. + # factor np.log(2) -> decay rate to half-life; factor 60 -> minutes to hours + hl_calc = ((np.log(2)/((0.16030 + 0.04018*(((inside_temp-273.15)-20.615)/10.585) +0.02176*(((humidity*100)-45.235)/28.665) -0.14369 -0.02636*((inside_temp-273.15)-20.615)/10.585)))/60) - if (hl_calc <= 0): - hl_calc = 6.43 - return hl_calc + + return np.where(hl_calc <= 0, 6.43, np.minimum(6.43, hl_calc)) Virus.types = { diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index b1959ca8..896c1c6b 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -179,12 +179,12 @@ def sr_model(): @pytest.mark.parametrize( ["exposed_time_interval", "expected_deposited_exposure"], [ - [(0., 1.), 49.60355486823376], - [(1., 1.01), 0.5876509666617122], - [(1.01, 1.02), 0.5726560233646302], - [(12., 12.01), 0.016288631412456556], - [(12., 24.), 716.6229828782525], - [(0., 24.), 777.8717785392312], + [(0., 1.), 42.63222033436878], + [(1., 1.01), 0.485377549596179], + [(1.01, 1.02), 0.47058239520823814], + [(12., 12.01), 0.01622776617499709], + [(12., 24.), 595.1115223695439], + [(0., 24.), 645.8401125684933], ] ) def test_exposure_model_integral_accuracy(exposed_time_interval, diff --git a/cara/tests/models/test_virus.py b/cara/tests/models/test_virus.py index 0f875fa8..83f4d1f0 100644 --- a/cara/tests/models/test_virus.py +++ b/cara/tests/models/test_virus.py @@ -7,18 +7,18 @@ from cara import models @pytest.mark.parametrize( "inside_temp, humidity, expected_halflife, expected_decay_constant", [ - [293.15, 0.5, 35.67710693238622, 0.01942834607844098], - [272.15, 0.7, 96.40459058258793, 0.007189981061805762], - [300.15, 1., 10.418034697539541, 0.0665333914393324], - [300.15, 0., 1.1, 0.6301338005090411], + [293.15, 0.5, 0.5947447349860315, 1.1654532436949188], + [272.15, 0.7, 1.6070844193207476, 0.4313072619127947], + [300.15, 1., 0.17367078830147223, 3.9911558376571805], + [300.15, 0., 6.43, 0.10779893943389507], [np.array([272.15, 300.15]), np.array([0.7, 0.]), - np.array([96.40459058258793, 1.1]), np.array([0.007189981061805762, 0.6301338005090411])], + np.array([1.60708442, 6.43]), np.array([0.43130726, 0.10779894])], [np.array([293.15, 300.15]), np.array([0.5, 1.]), - np.array([35.67710693238622, 10.418034697539541]), np.array([0.01942834607844098, 0.0665333914393324])] + np.array([0.59474473, 0.17367079]), np.array([1.16545324, 3.99115584])] ], ) def test_decay_constant(inside_temp, humidity, expected_halflife, expected_decay_constant): - npt.assert_equal(models.Virus.types['SARS_CoV_2'].halflife(humidity, inside_temp), - expected_halflife) - npt.assert_equal(models.Virus.types['SARS_CoV_2'].decay_constant(humidity, inside_temp), - expected_decay_constant) \ No newline at end of file + npt.assert_almost_equal(models.Virus.types['SARS_CoV_2'].halflife(humidity, inside_temp), + expected_halflife) + npt.assert_almost_equal(models.Virus.types['SARS_CoV_2'].decay_constant(humidity, inside_temp), + expected_decay_constant) \ No newline at end of file diff --git a/cara/tests/test_full_algorithm.py b/cara/tests/test_full_algorithm.py index 7afbd1ce..53a4e99e 100644 --- a/cara/tests/test_full_algorithm.py +++ b/cara/tests/test_full_algorithm.py @@ -84,12 +84,13 @@ class SimpleConcentrationModel: """ removal rate lambda in h^-1, excluding the deposition rate. """ + hl_calc = ((ln2/((0.16030 + 0.04018*(((293-273.15)-20.615)/10.585) + +0.02176*(((self.humidity*100)-45.235)/28.665) + -0.14369 + -0.02636*((293-273.15)-20.615)/10.585)))/60) return (self.lambda_ventilation - + ln2/(np.maximum(1.1, (0.693 / ((0.16030 + 0.04018 * (((21) - 20.615) / 10.585) - + 0.02176*(((self.humidity * 100) - 45.235) / 28.665) - - 0.14369 - - 0.02636*((21-20.615)/10.585))))))) + + ln2/(np.where(hl_calc <= 0, 6.43, np.minimum(6.43, hl_calc)))) @method_cache def deposition_removal_coefficient(self) -> float: diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 98a90acb..055c587b 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -43,7 +43,7 @@ def test_concentrations(baseline_concentration_model): concentrations = [baseline_concentration_model.concentration(float(t)) for t in ts] npt.assert_allclose( concentrations, - [0.000000e+00, 2.122276e+01, 1.240684e-12, 2.122276e+01, 7.253047e-26], + [0.000000e+00, 2.046096e+01, 3.846725e-13, 2.046096e+01, 7.231966e-27], rtol=1e-6 ) @@ -94,7 +94,7 @@ def test_r0(baseline_exposure_model): # expected r0 was computed with a trapezoidal integration, using # a mesh of 100'000 pts per exposed presence interval. r0 = baseline_exposure_model.reproduction_number() - npt.assert_allclose(r0, 783.490035) + npt.assert_allclose(r0, 771.380385) def test_periodic_window(baseline_periodic_window, baseline_room): @@ -381,8 +381,8 @@ def build_exposure_model(concentration_model, short_range_model): @pytest.mark.parametrize( "month, expected_deposited_exposure", [ - ['Jan', 401.300989], - ['Jun', 2420.383151], + ['Jan', 359.140499], + ['Jun', 1385.917562], ], ) def test_exposure_hourly_dep(month,expected_deposited_exposure, baseline_sr_model): @@ -402,8 +402,8 @@ def test_exposure_hourly_dep(month,expected_deposited_exposure, baseline_sr_mode @pytest.mark.parametrize( "month, expected_deposited_exposure", [ - ['Jan', 402.348745], - ['Jun', 2558.632473], + ['Jan', 359.983716], + ['Jun', 1439.267381], ], ) def test_exposure_hourly_dep_refined(month,expected_deposited_exposure, baseline_sr_model): diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index 7d75b2d3..3bce07ae 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -310,13 +310,13 @@ def waiting_room_mc(): @pytest.mark.parametrize( "mc_model, expected_pi, expected_new_cases, expected_dose, expected_ER", [ - ["shared_office_mc", 6.75, 0.20, 3.594, 809], - ["classroom_mc", 9.58, 1.82, 9.664, 5624], + ["shared_office_mc", 5.55, 0.17, 2.699, 809], + ["classroom_mc", 9.58, 1.82, 9.034, 5624], ["ski_cabin_mc", 16.0, 0.47, 17.315, 7966], - ["skagit_chorale_mc",72.20, 43.32, 134.234, 190422], - ["bus_ride_mc", 14.08, 9.43, 9.26, 5419], - ["gym_mc", 0.49, 0.14, 0.226, 1145], - ["waiting_room_mc", 2.02, 0.28, 1.104, 737], + ["skagit_chorale_mc",61.01, 36.53, 84.730, 190422], + ["bus_ride_mc", 10.59, 7.06, 6.65, 5419], + ["gym_mc", 0.43, 0.12, 0.197, 1145], + ["waiting_room_mc", 1.34, 0.18, 0.670, 737], ] ) def test_report_models(mc_model, expected_pi, expected_new_cases, @@ -337,10 +337,10 @@ def test_report_models(mc_model, expected_pi, expected_new_cases, @pytest.mark.parametrize( "mask_type, month, expected_pi, expected_dose, expected_ER", [ - ["No mask", "Jul", 11.81, 14.137, 809], - ["Type I", "Jul", 2.33, 1.305, 149], - ["FFP2", "Jul", 0.73, 0.351, 149], - ["Type I", "Feb", 0.62, 0.291, 149], + ["No mask", "Jul", 8.46, 8.113, 809], + ["Type I", "Jul", 1.44, 0.727, 149], + ["FFP2", "Jul", 0.43, 0.197, 149], + ["Type I", "Feb", 0.54, 0.253, 149], ], ) def test_small_shared_office_Geneva(mask_type, month, expected_pi, From 8b03425f00eb39cdd0a795d45b7ffc7769864127 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 9 May 2022 17:15:14 +0200 Subject: [PATCH 31/38] added retry decorator with default values --- cara/tests/test_monte_carlo_full_models.py | 2 ++ requirements.txt | 2 ++ 2 files changed, 4 insertions(+) diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index caed9c52..e5a25516 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -1,6 +1,7 @@ import numpy as np import numpy.testing as npt import pytest +from retry import retry import cara.monte_carlo as mc from cara import models,data @@ -309,6 +310,7 @@ def waiting_room_mc(): ) +@retry() @pytest.mark.parametrize( "mc_model, expected_pi, expected_new_cases, expected_dose, expected_ER", [ diff --git a/requirements.txt b/requirements.txt index 724bf5cc..a7a970cb 100644 --- a/requirements.txt +++ b/requirements.txt @@ -64,6 +64,7 @@ python-dateutil==2.8.2 pyzmq==22.1.0 requests==2.26.0 requests-unixsocket==0.2.0 +retry==0.9.2 scikit-learn==0.24.2 scipy==1.7.0 Send2Trash==1.7.1 @@ -76,6 +77,7 @@ threadpoolctl==2.2.0 timezonefinder==5.2.0 tornado==6.1 traitlets==5.0.5 +types-retry==0.9.7 urllib3==1.26.6 voila==0.2.10 wcwidth==0.2.5 From e3d9dc66ff50f4de57b440e67feaf2401344b2d2 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 9 May 2022 17:20:23 +0200 Subject: [PATCH 32/38] added retry requirement to setup.py --- setup.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/setup.py b/setup.py index e8136836..bc035ea7 100644 --- a/setup.py +++ b/setup.py @@ -30,10 +30,12 @@ REQUIREMENTS: dict = { 'numpy', 'psutil', 'python-dateutil', + 'retry', 'scipy', 'sklearn', 'timezonefinder', 'tornado', + 'types-retry', 'voila >=0.2.4', ], 'app': [], From 0b9582ce310a698338a631077ce1e47c8c8a46f6 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 10 May 2022 11:41:19 +0200 Subject: [PATCH 33/38] exposure of the inside_temp and humidity to the API --- cara/apps/calculator/model_generator.py | 15 +++++++++++---- cara/apps/calculator/static/js/form.js | 3 +++ cara/apps/templates/base/calculator.form.html.j2 | 2 ++ 3 files changed, 16 insertions(+), 4 deletions(-) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index 08cbc247..dd7f1e37 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -45,6 +45,7 @@ class FormData: floor_area: float hepa_amount: float hepa_option: bool + humidity: str infected_coffee_break_option: str #Used if infected_dont_have_breaks_with_exposed infected_coffee_duration: int #Used if infected_dont_have_breaks_with_exposed infected_dont_have_breaks_with_exposed: bool @@ -54,6 +55,7 @@ class FormData: infected_lunch_start: minutes_since_midnight #Used if infected_dont_have_breaks_with_exposed infected_people: int infected_start: minutes_since_midnight + inside_temp: float location_name: str location_latitude: float location_longitude: float @@ -100,6 +102,7 @@ class FormData: 'floor_area': 0., 'hepa_amount': 0., 'hepa_option': False, + 'humidity': '', 'infected_coffee_break_option': 'coffee_break_0', 'infected_coffee_duration': 5, 'infected_dont_have_breaks_with_exposed': False, @@ -109,6 +112,7 @@ class FormData: 'infected_lunch_start': '12:30', 'infected_people': _NO_DEFAULT, 'infected_start': '08:30', + 'inside_temp': 293., 'location_latitude': _NO_DEFAULT, 'location_longitude': _NO_DEFAULT, 'location_name': _NO_DEFAULT, @@ -240,11 +244,14 @@ class FormData: volume = self.room_volume else: volume = self.floor_area * self.ceiling_height - if self.room_heating_option: - humidity = 0.3 + if self.humidity == '': + if self.room_heating_option: + humidity = 0.3 + else: + humidity = 0.5 else: - humidity = 0.5 - room = models.Room(volume=volume, inside_temp=models.PiecewiseConstant((0, 24), (293,)), humidity=humidity) + humidity = float(self.humidity) + room = models.Room(volume=volume, inside_temp=models.PiecewiseConstant((0, 24), (self.inside_temp,)), humidity=humidity) infected_population = self.infected_population() diff --git a/cara/apps/calculator/static/js/form.js b/cara/apps/calculator/static/js/form.js index fd4f3c73..e4597acc 100644 --- a/cara/apps/calculator/static/js/form.js +++ b/cara/apps/calculator/static/js/form.js @@ -783,6 +783,9 @@ $(document).ready(function () { templateSelection: formatLocationSelection }); + // Logic for the API requests. Always set humity input as the empty string so that we can profit from the "room_heating_option default" values for humidity. + $("[name='humidity']").val(""); + function formatlocation(suggestedLocation) { // Function is called for each location from the geocoding API. diff --git a/cara/apps/templates/base/calculator.form.html.j2 b/cara/apps/templates/base/calculator.form.html.j2 index 499fc761..e1b66074 100644 --- a/cara/apps/templates/base/calculator.form.html.j2 +++ b/cara/apps/templates/base/calculator.form.html.j2 @@ -137,6 +137,8 @@ + +
From 96af7827370f342b11757992e33b0f4381f613a3 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 10 May 2022 11:56:49 +0200 Subject: [PATCH 34/38] Added humidity and inside_temp to the default values for the baseline form data --- cara/apps/calculator/model_generator.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index dd7f1e37..c7ebd410 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -693,6 +693,7 @@ def baseline_raw_form_data() -> typing.Dict[str, typing.Union[str, float]]: 'floor_area': '', 'hepa_amount': '250', 'hepa_option': '0', + 'humidity': '', 'infected_coffee_break_option': 'coffee_break_4', 'infected_coffee_duration': '10', 'infected_dont_have_breaks_with_exposed': '1', @@ -702,6 +703,7 @@ def baseline_raw_form_data() -> typing.Dict[str, typing.Union[str, float]]: 'infected_lunch_start': '12:30', 'infected_people': '1', 'infected_start': '09:00', + 'inside_temp': 293., 'location_latitude': 46.20833, 'location_longitude': 6.14275, 'location_name': 'Geneva', From 45a8e73a487f3b50431e7bece6057621aa355436 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 10 May 2022 17:02:40 +0200 Subject: [PATCH 35/38] added retry decorator on test_webapp.py --- cara/tests/apps/calculator/test_webapp.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/cara/tests/apps/calculator/test_webapp.py b/cara/tests/apps/calculator/test_webapp.py index e6b4eedf..fa5e37cc 100644 --- a/cara/tests/apps/calculator/test_webapp.py +++ b/cara/tests/apps/calculator/test_webapp.py @@ -2,11 +2,12 @@ from pathlib import Path import pytest import tornado.testing +from retry import retry import cara.apps.calculator from cara.apps.calculator.report_generator import generate_permalink -_TIMEOUT = 30. +_TIMEOUT = 20. @pytest.fixture def app(): @@ -42,6 +43,7 @@ async def test_404(http_server_client): assert resp.code == 404 +@retry() class TestBasicApp(tornado.testing.AsyncHTTPTestCase): def get_app(self): return cara.apps.calculator.make_app() @@ -70,6 +72,7 @@ class TestBasicApp(tornado.testing.AsyncHTTPTestCase): assert 'expected number of new cases is' in response.body.decode() +@retry() class TestCernApp(tornado.testing.AsyncHTTPTestCase): def get_app(self): cern_theme = Path(cara.apps.calculator.__file__).parent.parent / 'themes' / 'cern' @@ -82,6 +85,7 @@ class TestCernApp(tornado.testing.AsyncHTTPTestCase): assert 'expected number of new cases is' in response.body.decode() +retry() class TestOpenApp(tornado.testing.AsyncHTTPTestCase): def get_app(self): return cara.apps.calculator.make_app(calculator_prefix="/mycalc") From fe078f1446c3264cb34b21c54d76f0bd471351f3 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 10 May 2022 16:13:27 +0200 Subject: [PATCH 36/38] Added zenodo DOI for repository into the readme and index.html.j2 page --- README.md | 1 + cara/apps/templates/base/index.html.j2 | 1 + 2 files changed, 2 insertions(+) diff --git a/README.md b/README.md index e358a978..7d4d8778 100644 --- a/README.md +++ b/README.md @@ -39,6 +39,7 @@ Andre Henriques1, Luis Aleixo1, Marco Andreini1 **For the use of the CARA web app** CARA – COVID Airborne Risk Assessment tool +[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.6520432.svg)](https://doi.org/10.5281/zenodo.6520432) © Copyright 2020-2021 CERN. All rights not expressly granted are reserved. **For use of the model** diff --git a/cara/apps/templates/base/index.html.j2 b/cara/apps/templates/base/index.html.j2 index d6ee53d6..69a6511c 100644 --- a/cara/apps/templates/base/index.html.j2 +++ b/cara/apps/templates/base/index.html.j2 @@ -64,6 +64,7 @@ For use of the CARA web app:
  • CARA – COVID Airborne Risk Assessment tool
    • + DOI
    © Copyright 2020-2021 CERN. All rights not expressly granted are reserved.
    Licensed under the Apache License, Version 2.0
    LICENSE From 601ff171d62bdcb8d9f40f9814c32e68eaaa0e48 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 10 May 2022 16:16:43 +0200 Subject: [PATCH 37/38] Added newline in readme --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index 7d4d8778..12d3b59a 100644 --- a/README.md +++ b/README.md @@ -39,6 +39,7 @@ Andre Henriques1, Luis Aleixo1, Marco Andreini1 **For the use of the CARA web app** CARA – COVID Airborne Risk Assessment tool + [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.6520432.svg)](https://doi.org/10.5281/zenodo.6520432) © Copyright 2020-2021 CERN. All rights not expressly granted are reserved. From 26cf6fb0f0f1e50159d592fea10029c9bcdc670f Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 10 May 2022 17:19:59 +0200 Subject: [PATCH 38/38] added missing newline --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index 12d3b59a..b6354ce1 100644 --- a/README.md +++ b/README.md @@ -38,9 +38,11 @@ Andre Henriques1, Luis Aleixo1, Marco Andreini1 ### Reference and Citation **For the use of the CARA web app** + CARA – COVID Airborne Risk Assessment tool [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.6520432.svg)](https://doi.org/10.5281/zenodo.6520432) + © Copyright 2020-2021 CERN. All rights not expressly granted are reserved. **For use of the model**