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,