pingu_98/cara
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge branch 'model/humidity_function' into 'master'

Browse source 
Implement new virus half-life formula (humidity & temp. dependent)

See merge request cara/cara!351
This commit is contained in:
James Devine 2022-05-11 09:52:45 +02:00
commit 55c8c8d58c
15 changed files with 149 additions and 121 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

20
cara/apps/calculator/model_generator.py
View file

@ -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, 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()
@ -329,13 +336,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 +349,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,
@ -689,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',
@ -698,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',

3
cara/apps/calculator/static/js/form.js
View file

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

28
cara/apps/expert.py
View file

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

2
cara/apps/templates/base/calculator.form.html.j2
View file

@ -137,6 +137,8 @@
<input type="text" name="location_name" value="Geneva, CHE">
<input type="text" name="location_latitude" value="46.20833">
<input type="text" name="location_longitude" value="6.14275">
<input type="text" name="inside_temp" value="293">
<input type="text" name="humidity" value="">
</div>
<hr width="80%">

69
cara/models.py
View file

@ -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,9 @@ 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)
transitions.update(room.inside_temp.transition_times)
transitions.update(self.outside_temp.transition_times)
return transitions
@ -304,7 +304,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,28 +439,35 @@ 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:
# Viral inactivation per hour (h^-1) (function of humidity)
return np.log(2) / self.halflife(humidity)
def decay_constant(self, humidity: _VectorisedFloat, inside_temp: _VectorisedFloat) -> _VectorisedFloat:
# Viral inactivation per hour (h^-1) (function of humidity and inside temperature)
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. with correction https://doi.org/10.1080/02786826.2020.1829536
# 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)
return np.where(hl_calc <= 0, 6.43, np.minimum(6.43, hl_calc))
Virus.types = {
@ -917,9 +924,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.value(time))
+ self.ventilation.air_exchange(self.room, time)
)
@ -950,7 +957,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

7
cara/tests/apps/calculator/test_model_generator.py
View file

@ -60,7 +60,6 @@ def test_ventilation_slidingwindow(baseline_form: model_generator.FormData):
window = models.SlidingWindow(
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,
)
@ -92,7 +91,6 @@ def test_ventilation_hingedwindow(baseline_form: model_generator.FormData):
window = models.HingedWindow(
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,
)
@ -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=9),
inside_temp=models.PiecewiseConstant((0, 24), (293,)),
outside_temp=baseline_window.outside_temp,
window_height=1.6, opening_length=0.6,
)

3
cara/tests/conftest.py
View file

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

4
cara/tests/models/test_concentration_model.py
View file

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

14
cara/tests/models/test_exposure_model.py
View file

@ -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.), 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.), 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,

24
cara/tests/models/test_virus.py Normal file
View file

@ -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, 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([1.60708442, 6.43]), np.array([0.43130726, 0.10779894])],
[np.array([293.15, 300.15]), np.array([0.5, 1.]),
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_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)

11
cara/tests/test_full_algorithm.py
View file

@ -84,8 +84,13 @@ 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) )
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.where(hl_calc <= 0, 6.43, np.minimum(6.43, hl_calc))))
@method_cache
def deposition_removal_coefficient(self) -> float:
@ -461,7 +466,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,

40
cara/tests/test_known_quantities.py
View file

@ -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
@ -44,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.046096e+01, 3.846725e-13, 2.046096e+01, 7.231966e-27],
rtol=1e-6
)
@ -95,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, 771.380385)
def test_periodic_window(baseline_periodic_window, baseline_room):
@ -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,
@ -387,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', 359.140499],
['Jun', 1385.917562],
],
)
def test_exposure_hourly_dep(month,expected_deposited_exposure, baseline_sr_model):
@ -408,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', 359.983716],
['Jun', 1439.267381],
],
)
def test_exposure_hourly_dep_refined(month,expected_deposited_exposure, baseline_sr_model):

4
cara/tests/test_monte_carlo.py
View file

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

31
cara/tests/test_monte_carlo_full_models.py
View file

@ -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,
@ -312,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", 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],
["gym_mc", 0.45, 0.13, 0.208, 1145],
["waiting_room_mc", 1.59, 0.22, 0.821, 737],
["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",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,
@ -339,21 +337,20 @@ 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, 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,
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,
),

10
cara/tests/test_ventilation.py
View file

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