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

Adapted tests according to CARA paper

Browse source
This commit is contained in:
Luis Aleixo 2021-11-08 14:58:55 +01:00
parent 34ad2af425
commit 9f20172214
5 changed files with 183 additions and 115 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
cara/apps/expert.py
View file

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

39
cara/data/__init__.py
View file

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

2
cara/models.py
View file

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

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

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

243
cara/tests/test_monte_carlo_full_models.py
View file

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