pingu_98/cara
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Merge branch 'feature/integration_tests_with_mc' into 'feature/MonteCarlo'

Browse source 
Introducing integration tests for Monte-Carlo; new gravitational settlement

See merge request cara/cara!194
This commit is contained in:
Nicolas Mounet 2021-06-04 11:46:08 +00:00
commit 934ba106e8
6 changed files with 332 additions and 16 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
cara/models.py
View file

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

12
cara/tests/models/test_exposure_model.py
View file

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

12
cara/tests/test_known_quantities.py
View file

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

316
cara/tests/test_monte_carlo_full_models.py Normal file
View file

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

2
cara/tests/test_predefined_distributions.py
View file

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

2
cara/tests/test_sampleable_distribution.py
View file

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