Updated humidity value to be in percentage

2022-04-27 15:23:46 +02:00 · 2022-04-27 15:23:46 +02:00 · 5dcb698622
commit 5dcb698622
parent 71f36e0aa4
6 changed files with 50 additions and 26 deletions
--- 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))))
        
--- 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,
--- a/cara/tests/models/test_virus.py
+++ 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)
--- 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)))))))

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