updated half-life formula and respective tests

This commit is contained in:
Luis Aleixo 2022-05-09 16:29:52 +02:00
parent f518a893d5
commit 1e689af19c
6 changed files with 43 additions and 41 deletions

View file

@ -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 = {

View file

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

View file

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

View file

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

View file

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

View file

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