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

Merge branch 'feature/refined_expiration_model' into 'feature/MonteCarlo'

Browse source 
Implementation of BLO model for expiration

See merge request cara/cara!186
This commit is contained in:
Nicolas Mounet 2021-06-01 07:29:09 +00:00
commit 126ae3a846
7 changed files with 109 additions and 97 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

69
cara/models.py
View file

@ -37,6 +37,7 @@ import typing
import numpy as np
from scipy.interpolate import interp1d
import scipy.integrate
if not typing.TYPE_CHECKING:
from memoization import cached
@ -477,7 +478,7 @@ class Mask:
η_inhale: _VectorisedFloat
#: Global factor applied to filtration efficiency of masks when exhaling.
factor_exhale: _VectorisedFloat = 1.
factor_exhale: float = 1.
#: Pre-populated examples of Masks.
types: typing.ClassVar[typing.Dict[str, "Mask"]]
@ -529,31 +530,55 @@ class _ExpirationBase:
types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]]
def aerosols(self, mask: Mask):
# total volume of aerosols expired (cm^3).
"""
total volume of aerosols expired per volume of air (mL/cm^3).
"""
raise NotImplementedError("Subclass must implement")
@dataclass(frozen=True)
class Expiration(_ExpirationBase):
"""
Simple model based on four different sizes of particles emitted,
with different ejection factors. See Fig. 4 in L. Morawska et al,
Size distribution and sites of origin of droplets expelled from the
human respiratory tract during expiratory activities,
Aerosol Science 40 (2009) pp. 256 - 269.
BLO model for the expiration (G. Johnson et al., Modality of human
expired aerosol size distributions, Journal of Aerosol Science,
vol. 42, no. 12, pp. 839 – 851, 2011,
https://doi.org/10.1016/j.jaerosci.2011.07.009).
Here all diameters (d) are in microns.
"""
ejection_factor: typing.Tuple[float, ...]
particle_sizes: typing.Tuple[float, ...] = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4) # In cm.
#: factors assigned to resp. the B, L and O modes. They are
# charateristics of the kind of expiratory activity (e.g. breathing,
# speaking, singing, or shouting).
BLO_factors: typing.Tuple[float, float, float]
@cached()
def aerosols(self, mask: Mask):
def volume(diameter):
return (4 * np.pi * (diameter/2)**3) / 3
total = 0
for diameter, factor in zip(self.particle_sizes, self.ejection_factor):
contribution = (volume(diameter) * factor *
(1 - mask.exhale_efficiency(diameter*1e4)))
total += contribution
return total
""" Result is in mL.cm^-3 """
def volume(d):
return (np.pi * d**3) / 6.
def _Bmode(d: float) -> float:
# B-mode (see ref. above).
return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.262364)) *
np.exp(-1 * (np.log(d) - 0.989541) ** 2 / (2 * 0.262364 ** 2)))
def _Lmode(d: float) -> float:
# L-mode (see ref. above).
return ( (1 / d) * (1.0 / (np.sqrt(2 * np.pi) * 0.506818)) *
np.exp(-1 * (np.log(d) - 1.38629) ** 2 / (2 * 0.506818 ** 2)))
def _Omode(d: float) -> float:
# O-mode (see ref. above).
return ( (1 / d) * (0.0010008 / (np.sqrt(2 * np.pi) * 0.585005)) *
np.exp(-1 * (np.log(d) - 4.97673) ** 2 / (2 * 0.585005 ** 2)))
def integrand(d: float) -> float:
return (self.BLO_factors[0] * _Bmode(d) +
self.BLO_factors[1] * _Lmode(d) +
self.BLO_factors[2] * _Omode(d)
) * volume(d) * (1 - mask.exhale_efficiency(d))
# final result converted from microns^3/cm3 to mL/cm^3
return scipy.integrate.quad(integrand, 0.1, 30.)[0]*1e-12
@dataclass(frozen=True)
@ -581,11 +606,11 @@ class MultipleExpiration(_ExpirationBase):
_ExpirationBase.types = {
'Breathing': Expiration((0.084, 0.009, 0.003, 0.002)),
'Whispering': Expiration((0.11, 0.014, 0.004, 0.002)),
'Talking': Expiration((0.236, 0.068, 0.007, 0.011)),
'Unmodulated Vocalization': Expiration((0.751, 0.139, 0.0139, 0.059)),
'Superspreading event': Expiration((np.inf, np.inf, np.inf, np.inf)),
'Breathing': Expiration((1., 0., 0.)),
'Talking': Expiration((1., 1., 1.)),
'Shouting': Expiration((1., 5., 5.)),
'Singing': Expiration((1., 5., 5.)),
'Superspreading event': Expiration((np.inf, 0., 0.)),
}

12
cara/tests/conftest.py
View file

@ -7,11 +7,9 @@ import pytest
def baseline_model():
model = models.ConcentrationModel(
room=models.Room(volume=75),
ventilation=models.SlidingWindow(
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,
ventilation=models.AirChange(
active=models.SpecificInterval(((0,24),)),
air_exch=30.,
),
infected=models.InfectedPopulation(
number=1,
@ -19,7 +17,7 @@ def baseline_model():
presence=models.SpecificInterval(((0, 4), (5, 8))),
mask=models.Mask.types['No mask'],
activity=models.Activity.types['Light activity'],
expiration=models.Expiration.types['Unmodulated Vocalization'],
expiration=models.Expiration.types['Superspreading event'],
),
)
return model
@ -30,7 +28,7 @@ def baseline_exposure_model(baseline_model):
return models.ExposureModel(
baseline_model,
exposed=models.Population(
number=10,
number=1000,
presence=baseline_model.infected.presence,
activity=baseline_model.infected.activity,
mask=baseline_model.infected.mask,

8
cara/tests/models/test_concentration_model.py
View file

@ -14,7 +14,6 @@ from cara import models
{'air_change': np.array([100, 120])},
{'viral_load_in_sputum': np.array([5e8, 1e9])},
{'quantum_infectious_dose': np.array([50, 20])},
{'factor_exhale': np.array([0.92, 0.95])},
]
)
def test_concentration_model_vectorisation(override_params):
@ -24,7 +23,6 @@ def test_concentration_model_vectorisation(override_params):
'air_change': 100,
'viral_load_in_sputum': 1e9,
'quantum_infectious_dose': 50,
'factor_exhale': 0.95,
}
defaults.update(override_params)
@ -36,7 +34,7 @@ def test_concentration_model_vectorisation(override_params):
number=1,
presence=always,
mask=models.Mask(
factor_exhale=defaults['factor_exhale'],
factor_exhale=0.95,
η_inhale=0.3,
),
activity=models.Activity(
@ -47,9 +45,7 @@ def test_concentration_model_vectorisation(override_params):
viral_load_in_sputum=defaults['viral_load_in_sputum'],
quantum_infectious_dose=defaults['quantum_infectious_dose'],
),
expiration=models.Expiration(
ejection_factor=(0.084, 0.009, 0.003, 0.002),
),
expiration=models.Expiration((1., 0., 0.)),
)
)
concentrations = c_model.concentration(10)

30
cara/tests/test_expiration.py
View file

@ -9,7 +9,7 @@ from cara import models
def test_multiple_wrong_weight_size():
weights = (1., 2., 3.)
e_base = models.Expiration((0.084, 0.009, 0.003, 0.002))
e_base = models.Expiration((1., 0., 0.))
with pytest.raises(
ValueError,
match=re.escape("expirations and weigths should contain the"
@ -19,11 +19,25 @@ def test_multiple_wrong_weight_size():
def test_multiple():
weights = (1., 2.)
e1 = models.Expiration((0.03, 0.02, 0.01, 0.005))
e2 = models.Expiration((0.05, 0.04, 0.03, 0.01))
weights = (1., 1.)
e1 = models.Expiration((1., 0., 0.))
e2 = models.Expiration((4., 5., 5.))
e_expected = models.Expiration((2.5, 2.5, 2.5))
e = models.MultipleExpiration([e1, e2], weights)
assert e.aerosols(models.Mask.types['No mask']) == (
e1.aerosols(models.Mask.types['No mask'])/3. +
2*e2.aerosols(models.Mask.types['No mask'])/3.
)
mask = models.Mask.types['Type I']
npt.assert_almost_equal(e_expected.aerosols(mask), e.aerosols(mask))
# expected values obtained from another code
@pytest.mark.parametrize(
"BLO_weights, expected_aerosols",
[
[(1.,0.,0.), 1.38924e-12],
[(1.,1.,1.), 1.07129e-10],
[(1.,5.,5.), 5.30088e-10],
],
)
def test_expiration_aerosols(BLO_weights, expected_aerosols):
mask = models.Mask.types['No mask']
e = models.Expiration(BLO_weights)
npt.assert_allclose(e.aerosols(mask), expected_aerosols, rtol=1e-4)

8
cara/tests/test_infected_population.py
View file

@ -8,7 +8,6 @@ import cara.models
"override_params", [
{'viral_load_in_sputum': np.array([5e8, 1e9])},
{'quantum_infectious_dose': np.array([50, 20])},
{'factor_exhale': np.array([0.92, 0.95])},
{'exhalation_rate': np.array([0.75, 0.81])},
]
)
@ -16,7 +15,6 @@ def test_infected_population_vectorisation(override_params):
defaults = {
'viral_load_in_sputum': 1e9,
'quantum_infectious_dose': 50,
'factor_exhale': 0.95,
'exhalation_rate': 0.75,
}
defaults.update(override_params)
@ -26,7 +24,7 @@ def test_infected_population_vectorisation(override_params):
number=1,
presence=office_hours,
mask=cara.models.Mask(
factor_exhale=defaults['factor_exhale'],
factor_exhale=0.95,
η_inhale=0.3,
),
activity=cara.models.Activity(
@ -37,9 +35,7 @@ def test_infected_population_vectorisation(override_params):
viral_load_in_sputum=defaults['viral_load_in_sputum'],
quantum_infectious_dose=defaults['quantum_infectious_dose'],
),
expiration=cara.models.Expiration(
ejection_factor=(0.084, 0.009, 0.003, 0.002),
),
expiration=cara.models.Expiration((1., 0., 0.)),
)
emission_rate = infected.emission_rate(10)
assert isinstance(emission_rate, np.ndarray)

77
cara/tests/test_known_quantities.py
View file

@ -6,21 +6,12 @@ import cara.models as models
import cara.data as data
def test_no_mask_aerosols(baseline_model):
exp = models.Expiration.types['Unmodulated Vocalization']
npt.assert_allclose(
exp.aerosols(models.Mask.types['No mask']),
6.077541e-12,
rtol=1e-5,
)
def test_no_mask_emission_rate(baseline_model):
rate = 151.938514
def test_no_mask_superspeading_emission_rate(baseline_model):
expected_rate = 970.
npt.assert_allclose(
[baseline_model.infected.emission_rate(t) for t in [0, 1, 4, 4.5, 5, 8, 9]],
[0, rate, rate, 0, 0, rate, 0],
rtol=1e-5
[0, expected_rate, expected_rate, 0, 0, expected_rate, 0],
rtol=1e-12
)
@ -48,23 +39,24 @@ def baseline_periodic_hepa():
def test_concentrations(baseline_model):
# expected concentrations were computed analytically
ts = [0, 4, 5, 7, 10]
concentrations = [baseline_model.concentration(t) for t in ts]
npt.assert_allclose(
concentrations,
[0.000000e+00, 2.619538e-01, 1.146999e-04, 2.619537e-01, 5.022289e-08],
rtol=1e-5
[0.000000e+00, 0.4189594, 1.6422648e-14, 0.4189594, 6.4374587e-28],
rtol=1e-6
)
def test_smooth_concentrations(baseline_model):
# We don't care about the actual concentrations in this test, but rather
# that the curve itself is smooth.
dx = 0.1
dy_limit = dx * 2 # Anything more than this is a bit steep.
dx = 0.002
dy_limit = 0.2 # Anything more than this (in relative) is a bit steep.
ts = np.arange(0, 10, dx)
concentrations = [baseline_model.concentration(t) for t in ts]
assert np.abs(np.diff(concentrations)).max() < dy_limit
assert np.abs(np.diff(concentrations)).max()/np.mean(concentrations) < dy_limit
def build_model(interval_duration):
@ -80,7 +72,7 @@ def build_model(interval_duration):
presence=models.SpecificInterval(((0, 4), (5, 8))),
mask=models.Mask.types['No mask'],
activity=models.Activity.types['Light activity'],
expiration=models.Expiration.types['Unmodulated Vocalization'],
expiration=models.Expiration.types['Superspreading event'],
),
)
return model
@ -98,8 +90,8 @@ def test_concentrations_startup(baseline_model):
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.
p = baseline_exposure_model.infection_probability()
npt.assert_allclose(p, 89.001748)
r0 = baseline_exposure_model.reproduction_number()
npt.assert_allclose(r0, 973.535888)
def test_periodic_window(baseline_periodic_window, baseline_room):
@ -183,15 +175,6 @@ def test_multiple_ventilation_HEPA_HVAC_AirChange(volume, expected_value):
expected_value,rtol=1e-5)
def test_expiration_aerosols():
mask = models.Mask.types['Type I']
exp1 = models.Expiration((0.751, 0.139, 0.0139, 0.059),
particle_sizes = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4))
exp2 = models.Expiration((0.059, 0.0139, 0.751, 0.139),
particle_sizes = (5.5e-4, 3.5e-4, 0.8e-4, 1.8e-4))
npt.assert_allclose(exp1.aerosols(mask), exp2.aerosols(mask), rtol=1e-5)
@pytest.mark.parametrize(
"time, expected_value",
[
@ -239,7 +222,7 @@ def build_hourly_dependent_model(month, intervals_open=((7.5, 8.5),),
presence=models.SpecificInterval(intervals_presence_infected),
mask=models.Mask.types['No mask'],
activity=models.Activity.types['Light activity'],
expiration=models.Expiration.types['Unmodulated Vocalization'],
expiration=models.Expiration.types['Superspreading event'],
),
)
return model
@ -260,7 +243,7 @@ def build_constant_temp_model(outside_temp, intervals_open=((7.5, 8.5),)):
presence=models.SpecificInterval(((0, 4), (5, 7.5))),
mask=models.Mask.types['No mask'],
activity=models.Activity.types['Light activity'],
expiration=models.Expiration.types['Unmodulated Vocalization'],
expiration=models.Expiration.types['Superspreading event'],
),
)
return model
@ -287,7 +270,7 @@ def build_hourly_dependent_model_multipleventilation(month, intervals_open=((7.5
presence=models.SpecificInterval(((0, 4), (5, 7.5))),
mask=models.Mask.types['No mask'],
activity=models.Activity.types['Light activity'],
expiration=models.Expiration.types['Unmodulated Vocalization'],
expiration=models.Expiration.types['Superspreading event'],
),
)
return model
@ -370,16 +353,16 @@ def build_exposure_model(concentration_model):
)
# expected r0 were computed with a trapezoidal integration, using
# expected quanta were computed with a trapezoidal integration, using
# a mesh of 100'000 pts per exposed presence interval.
@pytest.mark.parametrize(
"month, expected_r0",
"month, expected_quanta",
[
['Jan', 86.258533],
['Jun', 99.972103],
['Jan', 10.136783],
['Jun', 41.800377],
],
)
def test_r0_hourly_dep(month,expected_r0):
def test_quanta_hourly_dep(month,expected_quanta):
m = build_exposure_model(
build_hourly_dependent_model(
month,
@ -387,19 +370,19 @@ def test_r0_hourly_dep(month,expected_r0):
intervals_presence_infected=((8, 12), (13, 17))
)
)
p = m.infection_probability()
npt.assert_allclose(p, expected_r0)
quanta = m.quanta_exposure()
npt.assert_allclose(quanta, expected_quanta)
# expected r0 were computed with a trapezoidal integration, using
# expected quanta were computed with a trapezoidal integration, using
# a mesh of 100'000 pts per exposed presence interval.
@pytest.mark.parametrize(
"month, expected_r0",
"month, expected_quanta",
[
['Jan', 86.422736],
['Jun', 99.982323],
['Jan', 10.19818],
['Jun', 44.130683],
],
)
def test_r0_hourly_dep_refined(month,expected_r0):
def test_quanta_hourly_dep_refined(month,expected_quanta):
m = build_exposure_model(
build_hourly_dependent_model(
month,
@ -408,5 +391,5 @@ def test_r0_hourly_dep_refined(month,expected_r0):
temperatures=data.GenevaTemperatures,
)
)
p = m.infection_probability()
npt.assert_allclose(p, expected_r0)
quanta = m.quanta_exposure()
npt.assert_allclose(quanta, expected_quanta)

2
cara/tests/test_monte_carlo.py
View file

@ -53,7 +53,7 @@ def baseline_mc_model() -> cara.monte_carlo.ConcentrationModel:
presence=cara.models.SpecificInterval(((0, 4), (5, 8))),
mask=cara.models.Mask.types['No mask'],
activity=cara.models.Activity.types['Light activity'],
expiration=cara.models.Expiration.types['Unmodulated Vocalization'],
expiration=cara.models.Expiration.types['Breathing'],
),
)
return mc_model