diff --git a/cara/models.py b/cara/models.py index 794ed21b..1561b09d 100644 --- a/cara/models.py +++ b/cara/models.py @@ -731,8 +731,8 @@ class ConcentrationModel: return self.infected.virus def infectious_virus_removal_rate(self, time: float) -> _VectorisedFloat: - # Particle deposition on the floor - vg = 1 * 10 ** -4 + # Particle deposition on the floor (value from CERN-OPEN-2021-04) + vg = 1.88e-4 # Height of the emission source to the floor - i.e. mouth/nose (m) h = 1.5 # Deposition rate (h^-1) diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index 35c5ac0c..1c6c6d3d 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -136,12 +136,12 @@ def conc_model(): # expected quanta were computed with a trapezoidal integration, using # a mesh of 10'000 pts per exposed presence interval. @pytest.mark.parametrize("exposed_time_interval, expected_quanta", [ - [(0, 1), 5.4869151], - [(1, 1.01), 0.064013521], - [(1.01, 1.02), 0.062266596], - [(12, 12.01), 0.0019025904], - [(12, 24), 78.190763], - [(0, 24), 84.866592], + [(0, 1), 5.3334352], + [(1, 1.01), 0.061759078], + [(1.01, 1.02), 0.060016487], + [(12, 12.01), 0.0019012647], + [(12, 24), 75.513005], + [(0, 24), 81.956988], ] ) def test_exposure_model_integral_accuracy(exposed_time_interval, diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 4af1c06f..ab24bb08 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -44,7 +44,7 @@ def test_concentrations(baseline_model): concentrations = [baseline_model.concentration(t) for t in ts] npt.assert_allclose( concentrations, - [0.000000e+00, 0.4189594, 1.6422648e-14, 0.4189594, 6.4374587e-28], + [0.000000e+00, 0.41611256, 1.3205628e-14, 0.41611256, 4.1909001e-28], rtol=1e-6 ) @@ -91,7 +91,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, 973.535888) + npt.assert_allclose(r0, 972.880852) def test_periodic_window(baseline_periodic_window, baseline_room): @@ -359,8 +359,8 @@ def build_exposure_model(concentration_model): @pytest.mark.parametrize( "month, expected_quanta", [ - ['Jan', 10.136783], - ['Jun', 41.800377], + ['Jan', 9.930854], + ['Jun', 37.962708], ], ) def test_quanta_hourly_dep(month,expected_quanta): @@ -379,8 +379,8 @@ def test_quanta_hourly_dep(month,expected_quanta): @pytest.mark.parametrize( "month, expected_quanta", [ - ['Jan', 10.19818], - ['Jun', 44.130683], + ['Jan', 9.989881], + ['Jun', 39.99636], ], ) def test_quanta_hourly_dep_refined(month,expected_quanta): diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py new file mode 100644 index 00000000..7c08757f --- /dev/null +++ b/cara/tests/test_monte_carlo_full_models.py @@ -0,0 +1,316 @@ +import numpy as np +import numpy.testing as npt +import pytest + +import cara.monte_carlo as mc +from cara import models,data +from cara.monte_carlo.data import activity_distributions, virus_distributions + +# TODO: seed better the random number generators +np.random.seed(2000) +SAMPLE_SIZE = 50000 +TOLERANCE = 0.05 + +# references values for infection_probability and expected new cases +# in the following tests, were obtained from the feature/mc branch + +@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) + """ + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=50, humidity=0.3), + ventilation=models.MultipleVentilation( + ( + 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, + ), + models.AirChange( + active=models.SpecificInterval(((0,24),)), + air_exch=0.25, + ), + ), + ), + infected=mc.InfectedPopulation( + number=1, + virus=virus_distributions['SARS_CoV_2_B117'], + presence=mc.SpecificInterval(((0, 2), (2.1, 4), (5, 7), (7.1, 9))), + mask=models.Mask(η_inhale=0.3), + activity=activity_distributions['Seated'], + expiration=models.MultipleExpiration( + expirations=(models.Expiration.types['Talking'], + models.Expiration.types['Breathing']), + weights=(0.3, 0.7)), + ), + ) + return mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=3, + presence=concentration_mc.infected.presence, + activity=models.Activity.types['Seated'], + mask=concentration_mc.infected.mask, + ), + ) + + +@pytest.fixture +def classroom_mc(): + """ + Corresponds to the 2nd line of Table 5 in CERN-OPEN-2021-04 + """ + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=160, humidity=0.3), + ventilation=models.MultipleVentilation( + ( + 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, + ), + models.AirChange( + active=models.SpecificInterval(((0,24),)), + air_exch=0.25, + ), + ), + ), + infected=mc.InfectedPopulation( + number=1, + virus=virus_distributions['SARS_CoV_2_B117'], + presence=mc.SpecificInterval(((0, 2), (2.5, 4), (5, 7), (7.5, 9))), + mask=models.Mask.types['No mask'], + activity=activity_distributions['Light activity'], + expiration=models.Expiration.types['Talking'], + ), + ) + return mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=19, + presence=concentration_mc.infected.presence, + activity=models.Activity.types['Seated'], + mask=concentration_mc.infected.mask, + ), + ) + + +@pytest.fixture +def ski_cabin_mc(): + """ + Corresponds to the 3rd line of Table 5 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, + ), + infected=mc.InfectedPopulation( + number=1, + virus=virus_distributions['SARS_CoV_2_B117'], + presence=mc.SpecificInterval(((0, 1/3),)), + mask=models.Mask(η_inhale=0.3), + activity=activity_distributions['Moderate activity'], + expiration=models.Expiration.types['Talking'], + ), + ) + return mc.ExposureModel( + 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, + ), + ) + + +@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. + """ + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=300, humidity=0.5), + ventilation=models.AirChange( + active=models.SpecificInterval(((0,24),)), + air_exch=6, + ), + infected=mc.InfectedPopulation( + number=2, + virus=virus_distributions['SARS_CoV_2_B117'], + presence=mc.SpecificInterval(((0, 1),)), + mask=models.Mask.types["No mask"], + activity=activity_distributions['Heavy exercise'], + expiration=models.Expiration.types['Breathing'], + ), + ) + return mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=28, + presence=concentration_mc.infected.presence, + activity=models.Activity.types['Heavy exercise'], + mask=concentration_mc.infected.mask, + ), + ) + + +@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%) + """ + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=100, humidity=0.5), + ventilation=models.AirChange( + active=models.SpecificInterval(((0,24),)), + air_exch=0.25, + ), + infected=mc.InfectedPopulation( + number=1, + virus=virus_distributions['SARS_CoV_2_B117'], + presence=mc.SpecificInterval(((0, 2),)), + mask=models.Mask.types["No mask"], + activity=activity_distributions['Seated'], + expiration=models.MultipleExpiration( + expirations=(models.Expiration.types['Talking'], + models.Expiration.types['Breathing']), + weights=(0.3, 0.7)), + ), + ) + return mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=14, + presence=concentration_mc.infected.presence, + activity=models.Activity.types['Seated'], + mask=concentration_mc.infected.mask, + ), + ) + + +@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=models.Expiration((5., 5., 5.)), + ), + ) + 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, + ), + ) + + +@pytest.mark.parametrize( + "mc_model, expected_pi, expected_new_cases, expected_dose, expected_qR", + [ + ["shared_office_mc", 10.7, 0.32, 0.954, 10.9], + ["classroom_mc", 36.1, 6.85, 13.0, 474.4], + ["ski_cabin_mc", 16.3, 0.49, 0.599, 123.4], + ["gym_mc", 2.25, 0.63, 0.01307, 16.4], + ["waiting_room_mc", 9.72, 1.36, 0.571, 58.9], + ["skagit_chorale_mc",29.9, 17.9, 1.90, 1414], + ] +) +def test_report_models(mc_model, expected_pi, expected_new_cases, + expected_dose, expected_qR, request): + mc_model = request.getfixturevalue(mc_model) + exposure_model = mc_model.build_model(size=SAMPLE_SIZE) + npt.assert_allclose(exposure_model.infection_probability().mean(), + expected_pi, rtol=TOLERANCE) + npt.assert_allclose(exposure_model.expected_new_cases().mean(), + expected_new_cases, rtol=TOLERANCE) + npt.assert_allclose(exposure_model.quanta_exposure().mean(), + expected_dose, rtol=TOLERANCE) + npt.assert_allclose( + exposure_model.concentration_model.infected.emission_rate_when_present().mean(), + expected_qR, rtol=TOLERANCE) + + +@pytest.mark.parametrize( + "mask_type, month, expected_pi, expected_dose, expected_qR", + [ + ["No mask", "Jul", 30.0, 6.764, 64.9], + ["Type I", "Jul", 10.2, 1.223, 11.7], + ["FFP2", "Jul", 4.0, 1.223, 11.7], + ["Type I", "Feb", 4.25, 0.357, 11.7], + ], +) +def test_small_shared_office_Geneva(mask_type, month, expected_pi, + expected_dose, expected_qR): + concentration_mc = mc.ConcentrationModel( + room=models.Room(volume=33, 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, + ), + models.AirChange( + active=models.SpecificInterval(((0,24),)), + air_exch=0.25, + ), + ), + ), + infected=mc.InfectedPopulation( + number=1, + virus=virus_distributions['SARS_CoV_2_B117'], + presence=mc.SpecificInterval(((9, 10+2/3), (10+5/6, 12.5), (13.5, 15+2/3), (15+5/6, 18))), + mask=models.Mask.types[mask_type], + activity=activity_distributions['Seated'], + expiration=models.MultipleExpiration( + expirations=(models.Expiration.types['Talking'], + models.Expiration.types['Breathing']), + weights=(0.33, 0.67)), + ), + ) + exposure_mc = mc.ExposureModel( + concentration_model=concentration_mc, + exposed=mc.Population( + number=1, + presence=concentration_mc.infected.presence, + activity=activity_distributions['Seated'], + mask=concentration_mc.infected.mask, + ), + ) + exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) + npt.assert_allclose(exposure_model.infection_probability().mean(), + expected_pi, rtol=TOLERANCE) + npt.assert_allclose(exposure_model.quanta_exposure().mean(), + expected_dose, rtol=TOLERANCE) + npt.assert_allclose( + exposure_model.concentration_model.infected.emission_rate_when_present().mean(), + expected_qR, rtol=TOLERANCE) diff --git a/cara/tests/test_predefined_distributions.py b/cara/tests/test_predefined_distributions.py index 1e066432..4197b2fb 100644 --- a/cara/tests/test_predefined_distributions.py +++ b/cara/tests/test_predefined_distributions.py @@ -5,7 +5,7 @@ import pytest from cara.monte_carlo.data import activity_distributions, virus_distributions # TODO: seed better the random number generators -np.random.seed(0) +np.random.seed(2000) # mean & std deviations from CERN-OPEN-2021-04 (Table 4) diff --git a/cara/tests/test_sampleable_distribution.py b/cara/tests/test_sampleable_distribution.py index 585b4dde..c30fcbf7 100644 --- a/cara/tests/test_sampleable_distribution.py +++ b/cara/tests/test_sampleable_distribution.py @@ -53,7 +53,7 @@ def test_lognormal(mean_gaussian, std_gaussian): assert len(samples) == sample_size npt.assert_allclose([samples.mean(), samples.std()], [exact_mean, exact_std], rtol=0.01) - npt.assert_allclose(selected_histogram, exact_dist, rtol=0.02) + npt.assert_allclose(selected_histogram, exact_dist, rtol=0.03) @pytest.mark.parametrize(