From 9f2017221459db6df6b1568f1640d0460b8829b3 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 8 Nov 2021 14:58:55 +0100 Subject: [PATCH] Adapted tests according to CARA paper --- cara/apps/expert.py | 2 +- cara/data/__init__.py | 39 +++- cara/models.py | 2 +- cara/tests/models/test_exposure_model.py | 12 +- cara/tests/test_monte_carlo_full_models.py | 243 ++++++++++++--------- 5 files changed, 183 insertions(+), 115 deletions(-) diff --git a/cara/apps/expert.py b/cara/apps/expert.py index a5f11f86..261afe75 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -498,7 +498,7 @@ baseline_model = models.ExposureModel( presence=models.SpecificInterval(((8., 12.), (13., 17.))), mask=models.Mask.types['No mask'], activity=models.Activity.types['Seated'], - expiration=models.Expiration.types['Talking'], + expiration=models.Expiration.types['Speaking'], host_immunity=0., ), evaporation_factor=0.3, diff --git a/cara/data/__init__.py b/cara/data/__init__.py index 353992e6..eea6f5c8 100644 --- a/cara/data/__init__.py +++ b/cara/data/__init__.py @@ -3,7 +3,7 @@ from cara import models # TODO: The values in this module to be removed and instead use the cara.data.weather functionality. -# average temperature of each month, hour per hour (from midnight to 11 pm) +# Geneva average temperature of each month, hour per hour (from midnight to 11 pm) Geneva_hourly_temperatures_celsius_per_hour = { 'Jan': [0.2, -0.3, -0.5, -0.9, -1.1, -1.4, -1.5, -1.5, -1.1, 0.1, 1.5, 2.8, 3.8, 4.4, 4.5, 4.4, 4.4, 3.9, 3.1, 2.7, 2.2, 1.7, 1.5, 1.1], @@ -31,6 +31,22 @@ Geneva_hourly_temperatures_celsius_per_hour = { 4.7, 5.2, 5.3, 5.2, 5.2, 4.7, 4.0, 3.7, 3.2, 2.8, 2.6, 2.2] } +# Toronto average temperature of each month, hour per hour (from midnight to 11 pm) +Toronto_hourly_temperatures_celsius_per_hour = { + "Jan": [ -2.9, -3.0, -3.2, -3.3, -3.3, -3.5, -3.7, -3.8, -3.9, -4.0, -4.1, -4.3, -4.3, -4.3, -4.1, -3.7, -3.2, -2.8, -2.6, -2.3, -2.2, -2.3, -2.6, -2.8], + "Feb": [ -2.4, -2.6, -2.8, -2.9, -2.9, -3.1, -3.3, -3.4, -3.6, -3.8, -3.9, -4.0, -4.3, -4.2, -3.7, -3.2, -2.6, -2.1, -1.7, -1.5, -1.3, -1.4, -1.6, -2.1], + "Mar": [ 1.3, 1.0, 0.7, 0.5, 0.4, 0.1, -0.0, -0.2, -0.4, -0.5, -0.7, -0.8, -0.9, -0.3, 0.4, 1.0, 1.6, 2.0, 2.3, 2.7, 2.7, 2.7, 2.4, 1.9], + "Apr": [ 6.8, 6.5, 6.3, 5.9, 5.7, 5.4, 5.1, 4.9, 4.6, 4.4, 4.2, 4.3, 4.8, 5.5, 6.1, 6.7, 7.1, 7.6, 7.8, 8.1, 8.2, 8.2, 8.0, 7.6 ], + "May": [ 13.0, 12.6, 12.2, 11.8, 11.5, 11.2, 10.8, 10.5, 10.2, 9.9, 9.8, 10.0, 10.9, 11.6, 12.2, 12.7, 13.2, 13.6, 13.9, 14.3, 14.4, 14.3, 14.2, 13.8], + "Jun": [ 18.9, 18.2, 17.8, 17.4, 17.0, 16.6, 16.2, 15.9, 15.6, 15.4, 15.3, 15.8, 16.5, 17.3, 17.9, 18.4, 18.9, 19.4, 19.7, 20.1, 20.3, 20.3, 20.1, 19.7], + "Jul": [ 22.1, 21.4, 20.9, 20.5, 20.0, 19.6, 19.1, 18.9, 18.6, 18.3, 18.1, 18.5, 19.4, 20.3, 21.0, 21.6, 22.2, 22.7, 23.1, 23.4, 23.6, 23.5, 23.3, 22.9], + "Aug": [ 22.0, 21.4, 21.0, 20.7, 20.3, 20.0, 19.6, 19.3, 19.1, 18.8, 18.5, 18.4, 19.4, 20.3, 21.1, 21.7, 22.2, 22.8, 23.1, 23.4, 23.5, 23.3, 23.1, 22.6], + "Sep": [ 18.2, 17.8, 17.4, 17.3, 17.0, 16.6, 16.3, 16.0, 15.8, 15.5, 15.4, 15.0, 15.6, 16.6, 17.5, 18.2, 18.7, 19.2, 19.6, 19.8, 19.7, 19.6, 19.2, 18.6], + "Oct": [ 11.1, 10.9, 10.6, 10.5, 10.2, 10.1, 9.8, 9.7, 9.5, 9.3, 9.2, 9.0, 9.0, 9.7, 10.5, 11.2, 11.7, 12.2, 12.4, 12.6, 12.6, 12.3, 11.8, 11.3], + "Nov": [ 5.3, 5.1, 5.0, 4.7, 4.6, 4.4, 4.3, 4.2, 4.1, 4.0, 3.9, 3.8, 3.7, 4.0, 4.6, 5.2, 5.7, 6.1, 6.2, 6.4, 6.3, 6.0, 5.5, 5.3], + "Dec": [ 0.4, 0.3, 0.2, 0.0, -0.1, -0.2, -0.4, -0.5, -0.6, -0.7, -0.8, -0.8, -0.9, -0.9, -0.6, -0.2, 0.3, 0.7, 0.9, 1.1, 1.1, 0.9, 0.6, 0.5] + } + # Geneva hourly temperatures as piecewise constant function (in Kelvin). GenevaTemperatures_hourly = { @@ -44,8 +60,27 @@ GenevaTemperatures_hourly = { } -# Same temperatures on a finer temperature mesh (every 6 minutes). +# Toronto hourly temperatures as piecewise constant function (in Kelvin). +TorontoTemperatures_hourly = { + month: models.PiecewiseConstant( + # NOTE: It is important that the time type is float, not np.float, in + # order to allow hashability (for caching). + tuple(float(time) for time in range(25)), + tuple(273.15 + np.array(temperatures)), + ) + for month, temperatures in Toronto_hourly_temperatures_celsius_per_hour.items() +} + + +# Same Geneva temperatures on a finer temperature mesh (every 6 minutes). GenevaTemperatures = { month: GenevaTemperatures_hourly[month].refine(refine_factor=10) for month, temperatures in Geneva_hourly_temperatures_celsius_per_hour.items() } + + +# Same Toronto temperatures on a finer temperature mesh (every 6 minutes). +TorontoTemperatures = { + month: TorontoTemperatures_hourly[month].refine(refine_factor=10) + for month, temperatures in Toronto_hourly_temperatures_celsius_per_hour.items() +} \ No newline at end of file diff --git a/cara/models.py b/cara/models.py index b7c1ff95..6d484be4 100644 --- a/cara/models.py +++ b/cara/models.py @@ -628,7 +628,7 @@ class MultipleExpiration(_ExpirationBase): # The correspondence with the BLO coefficients is given. _ExpirationBase.types = { 'Breathing': Expiration(1.3844), # corresponds to B/L/O coefficients of (1, 0, 0) - 'Talking': Expiration(5.8925), # corresponds to B/L/O coefficients of (1, 1, 1) + 'Speaking': Expiration(5.8925), # corresponds to B/L/O coefficients of (1, 1, 1) 'Shouting': Expiration(10.0411), # corresponds to B/L/O coefficients of (1, 5, 5) 'Singing': Expiration(10.0411), # corresponds to B/L/O coefficients of (1, 5, 5) } diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index bdde0129..5976e659 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -64,7 +64,7 @@ def known_concentrations(func): presence=halftime, mask=models.Mask.types['Type I'], activity=models.Activity.types['Standing'], - virus=models.Virus.types['SARS_CoV_2_B117'], + virus=models.Virus.types['SARS_CoV_2_ALPHA'], expiration=models.Expiration.types['Speaking'], host_immunity=0., ) @@ -76,19 +76,19 @@ def known_concentrations(func): @pytest.mark.parametrize( "population, cm, expected_exposure, expected_probability", [ [populations[1], known_concentrations(lambda t: 36.), - np.array([432, 432]), np.array([77.2191556943, 74.6803506895])], + np.array([432, 432]), np.array([67.9503762594, 65.2366759251])], [populations[2], known_concentrations(lambda t: 36.), - np.array([432, 432]), np.array([61.1470214407, 65.2366759251])], + np.array([432, 432]), np.array([51.6749232285, 55.6374622042])], [populations[0], known_concentrations(lambda t: np.array([36., 72.])), - np.array([432, 864]), np.array([65.2366759251, 87.9151129926])], + np.array([432, 864]), np.array([55.6374622042, 80.3196524031])], [populations[1], known_concentrations(lambda t: np.array([36., 72.])), - np.array([432, 864]), np.array([77.2191556943, 93.589153588])], + np.array([432, 864]), np.array([67.9503762594, 87.9151129926])], [populations[2], known_concentrations(lambda t: np.array([36., 72.])), - np.array([432, 864]), np.array([61.1470214407, 87.9151129926])], + np.array([432, 864]), np.array([51.6749232285, 80.3196524031])], ]) def test_exposure_model_ndarray(population, cm, expected_exposure, expected_probability): diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index 52dcf46b..7e83df25 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -4,8 +4,7 @@ import pytest import cara.monte_carlo as mc from cara import models,data -from cara.monte_carlo.data import activity_distributions, virus_distributions -from cara.monte_carlo.data import expiration_distribution, expiration_distributions +from cara.monte_carlo.data import activity_distributions, virus_distributions, expiration_distributions, infectious_dose_distribution, viable_to_RNA_ratio_distribution from cara.apps.calculator.model_generator import build_expiration # TODO: seed better the random number generators @@ -19,32 +18,29 @@ TOLERANCE = 0.05 @pytest.fixture def shared_office_mc(): """ - Corresponds to the 1st line of Table 5 in CERN-OPEN-2021-04, but - speaking 30% of the time (instead of 1/3) + Corresponds to the 1st line of Table 4 in CERN-OPEN-2021-04 """ concentration_mc = mc.ConcentrationModel( - room=models.Room(volume=50, humidity=0.3), - ventilation=models.MultipleVentilation( - ( + room=models.Room(volume=50, humidity=0.5), + ventilation = models.MultipleVentilation( + ventilations=( models.SlidingWindow( - active=models.PeriodicInterval(period=120, duration=10), - inside_temp=models.PiecewiseConstant((0., 24.), (293,)), - outside_temp=models.PiecewiseConstant((0., 24.), (283,)), - window_height=1.6, opening_length=0.6, + 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, ), - models.AirChange( - active=models.SpecificInterval(((0., 24.), )), - air_exch=0.25, - ), - ), + models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=0.25), + ) ), infected=mc.InfectedPopulation( number=1, - virus=virus_distributions['SARS_CoV_2_ALPHA'], - presence=mc.SpecificInterval(((0., 2.), (2.1, 4.), (5., 7.), (7.1, 9.))), - mask=models.Mask(η_inhale=0.3), + presence=mc.SpecificInterval(present_times = ((0, 3.5), (4.5, 9))), + virus=virus_distributions['SARS_CoV_2_DELTA'], + mask=models.Mask.types['No mask'], activity=activity_distributions['Seated'], - expiration=build_expiration({'Speaking': 0.3, 'Breathing': 0.7}), + expiration=build_expiration({'Speaking': 0.33, 'Breathing': 0.67}), host_immunity=0., ), evaporation_factor=0.3, @@ -53,42 +49,40 @@ def shared_office_mc(): concentration_model=concentration_mc, exposed=mc.Population( number=3, - presence=concentration_mc.infected.presence, - activity=models.Activity.types['Seated'], - mask=concentration_mc.infected.mask, + presence=mc.SpecificInterval(present_times = ((0, 3.5), (4.5, 9))), + activity=activity_distributions['Seated'], + mask=models.Mask.types['No mask'], host_immunity=0., - ), + ) ) @pytest.fixture def classroom_mc(): """ - Corresponds to the 2nd line of Table 5 in CERN-OPEN-2021-04 + Corresponds to the 2nd line of Table 4 in CERN-OPEN-2021-04 """ concentration_mc = mc.ConcentrationModel( room=models.Room(volume=160, humidity=0.3), - ventilation=models.MultipleVentilation( - ( + ventilation = models.MultipleVentilation( + ventilations=( models.SlidingWindow( - active=models.PeriodicInterval(period=120, duration=10), + active=models.PeriodicInterval(period=120, duration=120), inside_temp=models.PiecewiseConstant((0., 24.), (293,)), - outside_temp=models.PiecewiseConstant((0., 24.), (283,)), - window_height=1.6, opening_length=0.6, + outside_temp=data.TorontoTemperatures['Dec'], + window_height=1.6, + opening_length=0.2, ), - models.AirChange( - active=models.SpecificInterval(((0., 24.),)), - air_exch=0.25, - ), - ), + models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=0.25), + ) ), infected=mc.InfectedPopulation( number=1, + presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), virus=virus_distributions['SARS_CoV_2_ALPHA'], - presence=mc.SpecificInterval(((0., 2.), (2.5, 4.), (5., 7.), (7.5, 9.))), - mask=models.Mask.types['No mask'], + mask=models.Mask.types["No mask"], activity=activity_distributions['Light activity'], - expiration=expiration_distributions['Speaking'], + expiration=build_expiration('Speaking'), host_immunity=0., ), evaporation_factor=0.3, @@ -97,9 +91,9 @@ def classroom_mc(): concentration_model=concentration_mc, exposed=mc.Population( number=19, - presence=concentration_mc.infected.presence, - activity=models.Activity.types['Seated'], - mask=concentration_mc.infected.mask, + presence=models.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), + activity=activity_distributions['Seated'], + mask=models.Mask.types["No mask"], host_immunity=0., ), ) @@ -108,21 +102,20 @@ def classroom_mc(): @pytest.fixture def ski_cabin_mc(): """ - Corresponds to the 3rd line of Table 5 in CERN-OPEN-2021-04 + Corresponds to the 3rd line of Table 4 in CERN-OPEN-2021-04 """ concentration_mc = mc.ConcentrationModel( - room=models.Room(volume=10, humidity=0.5), - ventilation=models.AirChange( - active=models.SpecificInterval(((0., 24.),)), - air_exch=0, - ), + room=models.Room(volume=10, humidity=0.3), + ventilation=models.MultipleVentilation( + (models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=0.0), + models.AirChange(active=models.PeriodicInterval(period=120, duration=120), air_exch=0.25))), infected=mc.InfectedPopulation( number=1, - virus=virus_distributions['SARS_CoV_2_ALPHA'], - presence=mc.SpecificInterval(((0., 1/3),)), - mask=models.Mask(η_inhale=0.3), + presence=models.SpecificInterval(((0, 20/60),)), + virus=virus_distributions['SARS_CoV_2_DELTA'], + mask=models.Mask.types['No mask'], activity=activity_distributions['Moderate activity'], - expiration=expiration_distributions['Speaking'], + expiration=build_expiration('Speaking'), host_immunity=0., ), evaporation_factor=0.3, @@ -131,20 +124,96 @@ def ski_cabin_mc(): concentration_model=concentration_mc, exposed=mc.Population( number=3, - presence=concentration_mc.infected.presence, - activity=models.Activity.types['Moderate activity'], - mask=concentration_mc.infected.mask, + presence=models.SpecificInterval(((0, 20/60),)), + activity=activity_distributions['Moderate activity'], + mask=models.Mask.types['No mask'], host_immunity=0., ), ) +@pytest.fixture +def skagit_chorale_mc(): + """ + Corresponds to the 4th line of Table 4 in CERN-OPEN-2021-04, + assuming viral is 10**9 instead of a LogCustomKernel distribution. + """ + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=810, humidity=0.5), + ventilation=models.AirChange( + active=models.PeriodicInterval(period=120, duration=120), + air_exch=0.7), + infected=mc.InfectedPopulation( + number=1, + presence=models.SpecificInterval(((0, 2.5), )), + virus=mc.SARSCoV2( + viral_load_in_sputum=10**9, + infectious_dose=infectious_dose_distribution, + viable_to_RNA_ratio=viable_to_RNA_ratio_distribution, + transmissibility_factor=1., + ), + mask=models.Mask.types['No mask'], + activity=activity_distributions['Moderate activity'], + expiration=build_expiration('Shouting'), + host_immunity=0., + ), + evaporation_factor=0.3, + ) + return mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=60, + presence=models.SpecificInterval(((0, 2.5), )), + activity=activity_distributions['Moderate activity'], + mask=models.Mask.types['No mask'], + host_immunity=0., + ), + ) + + +@pytest.fixture +def bus_ride_mc(): + """ + Corresponds to the 5th line of Table 4 in CERN-OPEN-2021-04, + assuming viral is 5*10**8 instead of a LogCustomKernel distribution. + """ + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=45, humidity=0.5), + ventilation=models.AirChange( + active=models.PeriodicInterval(period=120, duration=120), + air_exch=1.25), + infected=mc.InfectedPopulation( + number=1, + presence=models.SpecificInterval(((0, 1.67), )), + virus=mc.SARSCoV2( + viral_load_in_sputum=5*10**8, + infectious_dose=infectious_dose_distribution, + viable_to_RNA_ratio=viable_to_RNA_ratio_distribution, + transmissibility_factor=1., + ), + mask=models.Mask.types['No mask'], + activity=activity_distributions['Seated'], + expiration=build_expiration('Speaking'), + host_immunity=0., + ), + evaporation_factor=0.3, + ) + return mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=67, + presence=models.SpecificInterval(((0, 1.67), )), + activity=activity_distributions['Seated'], + mask=models.Mask.types['No mask'], + host_immunity=0., + ), + ) + + @pytest.fixture def gym_mc(): """ - Corresponds to the 4th line of Table 5 in CERN-OPEN-2021-04, - but there the expected number of cases is wrongly reported as 0.56 - while it should be 0.63. + Gym model for testing """ concentration_mc = mc.ConcentrationModel( room=models.Room(volume=300, humidity=0.5), @@ -178,8 +247,7 @@ def gym_mc(): @pytest.fixture def waiting_room_mc(): """ - Corresponds to the 5th line of Table 5 in CERN-OPEN-2021-04, but - speaking 30% of the time (instead of 20%) + Waiting room model for testing """ concentration_mc = mc.ConcentrationModel( room=models.Room(volume=100, humidity=0.5), @@ -210,51 +278,16 @@ def waiting_room_mc(): ) -@pytest.fixture -def skagit_chorale_mc(): - """ - Corresponds to the 6th line of Table 5 in CERN-OPEN-2021-04, but - infection probability should be 29.8% instead of 32%, and number of - new cases 17.9 instead of 21. - """ - concentration_mc = mc.ConcentrationModel( - room=models.Room(volume=810, humidity=0.5), - ventilation=models.AirChange( - active=models.SpecificInterval(((0,24),)), - air_exch=0.7, - ), - infected=mc.InfectedPopulation( - number=1, - virus=virus_distributions['SARS_CoV_2'], - presence=mc.SpecificInterval(((0., 2.5),)), - mask=models.Mask.types["No mask"], - activity=activity_distributions['Light activity'], - expiration=expiration_distribution((5., 5., 5.)), - host_immunity=0., - ), - evaporation_factor=0.3, - ) - return mc.ExposureModel( - concentration_model=concentration_mc, - exposed=mc.Population( - number=60, - presence=concentration_mc.infected.presence, - activity=models.Activity.types['Moderate activity'], - mask=concentration_mc.infected.mask, - host_immunity=0., - ), - ) - - @pytest.mark.parametrize( "mc_model, expected_pi, expected_new_cases, expected_dose, expected_ER", [ - ["shared_office_mc", 2.3, 0.069, 12.18, 138], - ["classroom_mc", 16.1, 3.06, 151.28, 5907], - ["ski_cabin_mc", 4.4, 0.13, 7.23, 1507], - ["gym_mc", 0.57, 0.16, 0.4852, 1145], - ["waiting_room_mc", 2.02, 0.28, 7.23, 779], - ["skagit_chorale_mc",11.42, 6.85, 36.55, 28572], + ["shared_office_mc", 6.03, 0.18, 24.55, 809], + ["classroom_mc", 10.0, 2.0, 79.98, 5624], + ["ski_cabin_mc", 17.0, 0.5, 40.25, 7966], + ["skagit_chorale_mc",70, 42.5, 241.28, 190422], + ["bus_ride_mc", 12.0, 8.0, 63.79, 5419], + ["gym_mc", 0.45, 0.13, 0.4852, 1145], + ["waiting_room_mc", 1.59, 0.22, 7.23, 737], ] ) def test_report_models(mc_model, expected_pi, expected_new_cases, @@ -275,10 +308,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.68, 84.93, 838], - ["Type I", "Jul", 2.12, 15.41, 147], - ["FFP2", "Jul", 0.66, 15.62, 160], - ["Type I", "Feb", 0.73, 4.58, 151], + ["No mask", "Jul", 10.02, 84.54, 809], + ["Type I", "Jul", 1.7, 15.64, 149], + ["FFP2", "Jul", 0.51, 15.64, 149], + ["Type I", "Feb", 0.57, 4.59, 149], ], ) def test_small_shared_office_Geneva(mask_type, month, expected_pi,