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

Merge branch 'bugfix/short-range-f_inf' into 'master'

Browse source 
Fixed error with f_inf (short-range)

Closes #404

See merge request caimira/caimira!496
This commit is contained in:
Luis Aleixo 2024-07-15 09:16:33 +02:00
commit bad998c4c3
5 changed files with 220 additions and 73 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
caimira/apps/calculator/__init__.py
View file

@ -42,7 +42,7 @@ from .user import AuthenticatedUser, AnonymousUser
# calculator version. If the calculator needs to make breaking changes (e.g. change
# form attributes) then it can also increase its MAJOR version without needing to
# increase the overall CAiMIRA version (found at ``caimira.__version__``).
__version__ = "4.16.0"
__version__ = "4.16.1"
LOG = logging.getLogger("Calculator")

5
caimira/apps/calculator/report_generator.py
View file

@ -279,11 +279,14 @@ def uncertainties_plot(infection_probability: models._VectorisedFloat,
lower_percentiles: models._VectorisedFloat,
upper_percentiles: models._VectorisedFloat):
fig, axs = plt.subplots(2, 3,
fig, axes = plt.subplots(2, 3,
gridspec_kw={'width_ratios': [5, 0.5] + [1],
'height_ratios': [3, 1], 'wspace': 0},
sharey='row',
sharex='col')
# Type hint for axs
axs: np.ndarray = np.array(axes)
for y, x in [(0, 1)] + [(1, i + 1) for i in range(2)]:
axs[y, x].axis('off')

140
caimira/models.py
View file

@ -687,13 +687,8 @@ class _ExpirationBase:
def aerosols(self, mask: Mask):
"""
Total volume of aerosols expired per volume of exhaled air (mL/cm^3).
"""
raise NotImplementedError("Subclass must implement")
def jet_origin_concentration(self):
"""
Concentration of viruses at the jet origin (mL/m3).
Total volume of aerosols expired per volume of exhaled air
considering the outward mask efficiency (mL/cm^3).
"""
raise NotImplementedError("Subclass must implement")
@ -723,8 +718,9 @@ class Expiration(_ExpirationBase):
@cached()
def aerosols(self, mask: Mask):
"""
Total volume of aerosols expired per volume of exhaled air.
Result is in mL.cm^-3
Total volume of aerosols expired per volume
of exhaled air considering the outward mask
efficiency. Result is in mL.cm^-3.
"""
def volume(d):
return (np.pi * d**3) / 6.
@ -733,14 +729,6 @@ class Expiration(_ExpirationBase):
return self.cn * (volume(self.diameter) *
(1 - mask.exhale_efficiency(self.diameter))) * 1e-12
@cached()
def jet_origin_concentration(self):
def volume(d):
return (np.pi * d**3) / 6.
# Final result converted from microns^3/cm3 to mL/m3
return self.cn * volume(self.diameter) * 1e-6
@dataclass(frozen=True)
class MultipleExpiration(_ExpirationBase):
@ -889,9 +877,9 @@ class _PopulationWithVirus(Population):
def emission_rate_per_aerosol_per_person_when_present(self) -> _VectorisedFloat:
"""
The emission rate of virions in the expired air per mL of respiratory fluid,
per person, if the infected population is present, in (virion.cm^3)/(mL.h).
This method includes only the diameter-independent variables within the emission rate.
The emission rate of infectious respiratory particles (IRP) in the expired air per
mL of respiratory fluid, if the infected population is present, in (virions.cm^3)/(mL.h).
This method returns only the diameter-independent variables within the emission rate.
It should not be a function of time.
"""
raise NotImplementedError("Subclass must implement")
@ -900,12 +888,12 @@ class _PopulationWithVirus(Population):
def emission_rate_per_person_when_present(self) -> _VectorisedFloat:
"""
The emission rate if the infected population is present, per person
(in virions / h).
(in virions/h).
"""
return (self.emission_rate_per_aerosol_per_person_when_present() *
self.aerosols())
def emission_rate(self, time) -> _VectorisedFloat:
def emission_rate(self, time: float) -> _VectorisedFloat:
"""
The emission rate of the population vs time.
"""
@ -945,13 +933,13 @@ class EmittingPopulation(_PopulationWithVirus):
@method_cache
def emission_rate_per_aerosol_per_person_when_present(self) -> _VectorisedFloat:
"""
The emission rate of virions in the expired air per mL of respiratory fluid,
per person, if the infected population is present, in (virion.cm^3)/(mL.h).
This method includes only the diameter-independent variables within the emission rate.
The emission rate of infectious respiratory particles (IRP) in the expired air per
mL of respiratory fluid, if the infected population is present, in (virions.cm^3)/(mL.h).
This method returns only the diameter-independent variables within the emission rate.
It should not be a function of time.
"""
return self.known_individual_emission_rate
@dataclass(frozen=True)
class InfectedPopulation(_PopulationWithVirus):
@ -974,17 +962,19 @@ class InfectedPopulation(_PopulationWithVirus):
@method_cache
def emission_rate_per_aerosol_per_person_when_present(self) -> _VectorisedFloat:
"""
The emission rate of virions in the expired air per mL of respiratory fluid,
if the infected population is present, in (virion.cm^3)/(mL.h).
This method includes only the diameter-independent variables within the emission rate.
The emission rate of infectious respiratory particles (IRP) in the expired air per
mL of respiratory fluid, if the infected population is present, in (virions.cm^3)/(mL.h).
This method returns only the diameter-independent variables within the emission rate.
It should not be a function of time.
"""
# Note on units: exhalation rate is in m^3/h -> 1e6 conversion factor
# Returns the emission rate times the number of infected hosts in the room
# Conversion factor explanation:
# The exhalation rate is in m^3/h, therefore the 1e6 conversion factor
# is to convert m^3/h into cm^3/h to return (virions.cm^3)/(mL.h),
# so that we can then multiply by aerosols (mL/cm^3).
ER = (self.virus.viral_load_in_sputum *
self.activity.exhalation_rate *
self.fraction_of_infectious_virus() *
10 ** 6)
10**6)
return ER
@property
@ -1386,29 +1376,38 @@ class ShortRangeModel:
xstar[distances >= xstar])/𝛽r1/(xstar[distances >= xstar]
+ x0))**3
return factors
def _normed_jet_origin_concentration(self) -> _VectorisedFloat:
"""
The initial jet concentration at the source origin (mouth/nose), normalized by
normalization_factor in the ShortRange class (corresponding to the diameter-independent
variables). Results in mL.cm^-3.
"""
# The short range origin concentration does not consider the mask contribution.
return self.expiration.aerosols(mask=Mask.types['No mask'])
def _long_range_normed_concentration(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat:
"""
Virus long-range exposure concentration normalized by the
virus viral load, as function of time.
Virus long-range exposure concentration normalized by normalization_factor in the
ShortRange class, as function of time. Results in mL.cm^-3.
"""
return (concentration_model.concentration(time) /
concentration_model.virus.viral_load_in_sputum)
return (concentration_model.concentration(time) / self.normalization_factor(concentration_model.infected))
def _normed_concentration(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat:
"""
Virus short-range exposure concentration, as a function of time.
If the given time falls within a short-range interval it returns the
short-range concentration normalized by the virus viral load. Otherwise
it returns 0.
short-range concentration normalized by normalization_factor in the
Short-range class. Otherwise it returns 0. Results in mL.cm^-3.
"""
start, stop = self.presence.boundaries()[0]
# Verifies if the given time falls within a short-range interaction
if start <= time <= stop:
dilution = self.dilution_factor()
jet_origin_concentration = self.expiration.jet_origin_concentration()
# Long-range concentration normalized by the virus viral load
# Jet origin concentration normalized by the emission rate (except the BR)
normed_jet_origin_concentration = self._normed_jet_origin_concentration()
# Long-range concentration normalized by the emission rate (except the BR)
long_range_normed_concentration = self._long_range_normed_concentration(concentration_model, time)
# The long-range concentration values are then approximated using interpolation:
@ -1420,15 +1419,32 @@ class ShortRangeModel:
# Short-range concentration formula. The long-range concentration is added in the concentration method (ExposureModel).
# based on continuum model proposed by Jia et al (2022) - https://doi.org/10.1016/j.buildenv.2022.109166
return ((1/dilution)*(jet_origin_concentration - long_range_normed_concentration_interpolated))
return ((1/dilution)*(normed_jet_origin_concentration - long_range_normed_concentration_interpolated))
return 0.
def normalization_factor(self, infected: InfectedPopulation) -> _VectorisedFloat:
"""
The normalization factor applied to the short-range results. It refers to the emission
rate per aerosol without accounting for the exhalation rate (viral load and f_inf).
Result in (virions.cm^3)/(mL.m^3).
"""
# Re-use the emission rate method divided by the BR contribution.
return infected.emission_rate_per_aerosol_per_person_when_present() / infected.activity.exhalation_rate
def jet_origin_concentration(self, infected: InfectedPopulation) -> _VectorisedFloat:
"""
The initial jet concentration at the source origin (mouth/nose).
Returns the full result with the diameter dependent and independent variables, in virions/m^3.
"""
return self._normed_jet_origin_concentration() * self.normalization_factor(infected)
def short_range_concentration(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat:
"""
Virus short-range exposure concentration, as a function of time.
Factor of normalization applied back here. Results in virions/m^3.
"""
return (self._normed_concentration(concentration_model, time) *
concentration_model.virus.viral_load_in_sputum)
return (self._normed_concentration(concentration_model, time) *
self.normalization_factor(concentration_model.infected))
@method_cache
def _normed_short_range_concentration_cached(self, concentration_model: ConcentrationModel, time: float) -> _VectorisedFloat:
@ -1462,16 +1478,16 @@ class ShortRangeModel:
return start, stop
def _normed_jet_exposure_between_bounds(self,
concentration_model: ConcentrationModel,
time1: float, time2: float):
"""
Get the part of the integrated short-range concentration of
viruses in the air, between the times start and stop, coming
from the jet concentration, normalized by the viral load, and
without dilution.
from the jet concentration, normalized by normalization_factor,
and without dilution.
"""
start, stop = self.extract_between_bounds(time1, time2)
jet_origin = self.expiration.jet_origin_concentration()
# Note the conversion factor mL.cm^-3 -> mL.m^-3
jet_origin = self._normed_jet_origin_concentration() * 10**6
return jet_origin * (stop - start)
def _normed_interpolated_longrange_exposure_between_bounds(
@ -1479,20 +1495,24 @@ class ShortRangeModel:
time1: float, time2: float):
"""
Get the part of the integrated short-range concentration due
to the background concentration, normalized by the viral load
and the breathing rate, and without dilution.
to the background concentration, normalized by normalization_factor
together with breathing rate, and without dilution.
One needs to interpolate the integrated long-range concentration
for the particle diameters defined here.
TODO: make sure any potential extrapolation has a
negligible effect.
"""
start, stop = self.extract_between_bounds(time1, time2)
if stop<=start:
return 0.
# Note that for the correct integration one needs to isolate those parameters
# that are diameter-dependent from those that are diameter independent.
# Therefore, the diameter-independent parameters (viral load, f_ind and BR)
# are removed for the interpolation, and added back once the integration over
# the new aerosol diameters (done with the mean) is completed.
normed_int_concentration = (
concentration_model.integrated_concentration(start, stop)
/concentration_model.virus.viral_load_in_sputum
/concentration_model.infected.fraction_of_infectious_virus()
/concentration_model.infected.activity.exhalation_rate
)
normed_int_concentration_interpolated = np.interp(
@ -1566,7 +1586,8 @@ class CO2DataModel:
@dataclass(frozen=True)
class ExposureModel:
"""
Represents the exposure to a concentration of virions in the air.
Represents the exposure to a concentration of
infectious respiratory particles (IRP) in the air.
"""
data_registry: DataRegistry
@ -1693,7 +1714,6 @@ class ExposureModel:
self.exposed.activity.inhalation_rate *
(1 - self.exposed.mask.inhale_efficiency()))
# In the end we multiply the final results by the fraction of infectious virus of the vD equation.
return deposited_exposure
def deposited_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat:
@ -1709,8 +1729,7 @@ class ExposureModel:
deposited_exposure: _VectorisedFloat = 0.
for interaction in self.short_range:
start, stop = interaction.extract_between_bounds(time1, time2)
short_range_jet_exposure = interaction._normed_jet_exposure_between_bounds(
self.concentration_model, start, stop)
short_range_jet_exposure = interaction._normed_jet_exposure_between_bounds(start, stop)
short_range_lr_exposure = interaction._normed_interpolated_longrange_exposure_between_bounds(
self.concentration_model, start, stop)
dilution = interaction.dilution_factor()
@ -1741,11 +1760,12 @@ class ExposureModel:
interaction.activity.inhalation_rate
/dilution)
# Then we multiply by diameter-independent quantities: viral load
# and fraction of infected virions
# Then we multiply by the emission rate without the BR contribution (and conversion factor),
# and parameters of the vD equation (i.e. n_in).
deposited_exposure *= (
self.concentration_model.virus.viral_load_in_sputum
* (1 - self.exposed.mask.inhale_efficiency()))
(self.concentration_model.infected.emission_rate_per_aerosol_per_person_when_present() / (
self.concentration_model.infected.activity.exhalation_rate * 10**6)) *
(1 - self.exposed.mask.inhale_efficiency()))
# Long-range concentration
deposited_exposure += self.long_range_deposited_exposure_between_bounds(time1, time2)

8
caimira/tests/models/test_short_range_model.py
View file

@ -49,7 +49,7 @@ def test_short_range_model_ndarray(concentration_model, short_range_model):
model = short_range_model.build_model(SAMPLE_SIZE)
assert isinstance(model._normed_concentration(concentration_model, 10.75), np.ndarray)
assert isinstance(model.short_range_concentration(concentration_model, 10.75), np.ndarray)
assert isinstance(model._normed_jet_exposure_between_bounds(concentration_model, 10.75, 10.85), np.ndarray)
assert isinstance(model._normed_jet_exposure_between_bounds(10.75, 10.85), np.ndarray)
assert isinstance(model._normed_interpolated_longrange_exposure_between_bounds(concentration_model, 10.75, 10.85), np.ndarray)
assert isinstance(model.short_range_concentration(concentration_model, 14.0), float)
@ -106,9 +106,9 @@ def test_extract_between_bounds(short_range_model, time1, time2,
@pytest.mark.parametrize(
"time, expected_short_range_concentration", [
[8.5, 0.],
[10.5, 11.266605],
[10.6, 11.266605],
[11.0, 11.266605],
[10.5, 5.6333025],
[10.6, 5.6333025],
[11.0, 5.6333025],
[12.0, 0.],
]
)

138
caimira/tests/test_full_algorithm.py
View file

@ -16,6 +16,7 @@ from caimira.monte_carlo.data import (expiration_distributions,
expiration_BLO_factors,short_range_expiration_distributions,
short_range_distances,virus_distributions,activity_distributions)
SAMPLE_SIZE = 1_000_000
TOLERANCE = 0.04
@ -55,7 +56,7 @@ class SimpleConcentrationModel:
#: Number of infected people
num_infected: int = 1
#: Fraction of infected viruses (viable to RNA ratio)
#: Viable to RNA ratio
viable_to_RNA: _VectorisedFloat = 0.5
#: Host immunity factor (0. for not immune)
@ -97,6 +98,12 @@ class SimpleConcentrationModel:
return (self.lambda_ventilation
+ ln2/(np.where(hl_calc <= 0, 6.43, np.minimum(6.43, hl_calc))))
def fraction_of_infectious_virus(self) -> _VectorisedFloat:
"""
The fraction of infectious virus.
"""
return self.viable_to_RNA * (1 - self.HI)
@method_cache
def deposition_removal_coefficient(self) -> float:
"""
@ -181,8 +188,8 @@ class SimpleConcentrationModel:
return ( ( (0 if not self.infected_presence.triggered(t)
else self.f(lambda_rate,0))
+ result * np.exp(-lambda_rate*(t-ti)) )
* self.num_infected * self.viable_to_RNA
* (1. - self.HI) / self.room_volume)
* self.num_infected * self.fraction_of_infectious_virus()
/ self.room_volume)
@dataclass(frozen=True)
@ -263,6 +270,7 @@ class SimpleShortRangeModel:
we perform the integral of Np(d)*V(d) over diameter analytically
"""
vl = conc_model.viral_load
viable_to_RNA = conc_model.viable_to_RNA
dmin = self.diameter_min
dmax = self.diameter_max
result = 0.
@ -273,7 +281,7 @@ class SimpleShortRangeModel:
ymax = (np.log(dmax)-mu)/(sqrt2*sigma)-3.*sigma/sqrt2
result += ( (cn * famp * d0**3)/2. * np.exp(9*sigma**2/2.) *
(erf(ymax) - erf(ymin)) )
return vl * 1e-6 * result * np.pi/6.
return vl * viable_to_RNA * 1e-6 * result * np.pi/6.
def concentration(self, conc_model: SimpleConcentrationModel, time: float) -> _VectorisedFloat:
"""
@ -410,8 +418,8 @@ class SimpleExposureModel(SimpleConcentrationModel):
else self.f_with_fdep(lambda_rate,0,evaporation)*(t2-t1))
+ (primitive(t2) * np.exp(-lambda_rate*(t2-ti)) -
primitive(t1) * np.exp(-lambda_rate*(t1-ti)) ) )
* self.num_infected * self.viable_to_RNA
* (1. - self.HI) / self.room_volume)
* self.num_infected * self.fraction_of_infectious_virus()
/ self.room_volume)
@method_cache
def integrated_shortrange_concentration(self) -> _VectorisedFloat:
@ -430,7 +438,7 @@ class SimpleExposureModel(SimpleConcentrationModel):
res = (quad(integrand,
sr_model.diameter_min,sr_model.diameter_max,
epsabs=0.,limit=500)[0]
* self.viral_load * 1e-6 * (t2-t1) )
* self.viral_load * self.fraction_of_infectious_virus() * 1e-6 * (t2-t1) )
result += sr_model.breathing_rate * (
res-self.integrated_longrange_concentration(t1,t2,evaporation)
)/sr_model.dilution_factor()
@ -842,3 +850,119 @@ def test_exposure_with_shortrange_and_distributions(expo_sr_model_distr,
rtol=0.03
)
def exposure_model_from_parameter(data_registry, short_range_models, f_inf=0.5, viral_load=1e9, BR=1.25):
virus: models.SARSCoV2 = models.SARSCoV2(
viral_load_in_sputum=viral_load,
infectious_dose=50,
viable_to_RNA_ratio=f_inf,
transmissibility_factor=0.51,
)
c_model = mc.ConcentrationModel(
data_registry=data_registry,
room=models.Room(volume=50, humidity=0.3),
ventilation=models.AirChange(active=models.PeriodicInterval(period=120, duration=120),
air_exch=10_000_000),
infected=mc.InfectedPopulation(
data_registry=data_registry,
number=1,
presence=models.SpecificInterval(present_times=((8.5, 12.5), (13.5, 17.5))),
virus=virus,
mask=models.Mask.types['No mask'],
activity=models.Activity.types['Seated'],
expiration=expiration_distributions(data_registry)['Breathing'],
host_immunity=0.,
),
evaporation_factor=0.3,
)
return mc.ExposureModel(
data_registry=data_registry,
concentration_model=c_model,
short_range=short_range_models,
exposed=mc.Population(
number=1,
presence=models.SpecificInterval(present_times=((8.5, 12.5), (13.5, 17.5))),
mask=models.Mask.types['No mask'],
activity=models.Activity(inhalation_rate=BR, exhalation_rate=1.25),
host_immunity=0.,
),
geographical_data=models.Cases(),
).build_model(SAMPLE_SIZE)
@retry(tries=10)
def test_exposure_scale_with_f_inf(data_registry, sr_models):
"""
Exposure scaling test for the fraction of infectious virus.
"""
e_model_1: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=sr_models, f_inf=0.5)
e_model_2: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=sr_models, f_inf=1)
np.testing.assert_allclose(
2*e_model_1.deposited_exposure().mean(),
e_model_2.deposited_exposure().mean(), rtol=0.02
)
@retry(tries=10)
def test_exposure_scale_with_viral_load(data_registry, sr_models):
"""
Exposure scaling test for the viral load.
"""
e_model_1: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=sr_models, viral_load=1e9)
e_model_2: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=sr_models, viral_load=2e9)
np.testing.assert_allclose(
2*e_model_1.deposited_exposure().mean(),
e_model_2.deposited_exposure().mean(), rtol=0.02
)
@retry(tries=10)
def test_lr_exposure_scale_with_breathing_rate(data_registry):
"""
Exposure scaling test for the breathing rate when there are only long-range
interactions defined. Only the inhalation rate of the infected takes place
at the deposited exposure level.
"""
e_model_1: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=(), BR=1.25)
e_model_2: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=(), BR=2.5)
np.testing.assert_allclose(
2*e_model_1.deposited_exposure().mean(),
e_model_2.deposited_exposure().mean(), rtol=0.02
)
@retry(tries=10)
def test_exposure_scale_with_breathing_rate(data_registry, sr_models):
"""
Exposure scaling test for the breathing rate when long- and short-range
interactions are defined. We need to apply the multiplication factor
to the inhalation rate of the infected (long-range), but also for
each short-range interaction.
"""
e_model_1: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=sr_models, BR=1.25)
seated_act = models.Activity.types['Seated']
heavy_exercise_act = models.Activity.types['Heavy exercise']
sr_models_activity = (
mc.ShortRangeModel(
data_registry = data_registry,
expiration = short_range_expiration_distributions(data_registry)['Speaking'],
activity = models.Activity(inhalation_rate=seated_act.inhalation_rate * 2,
exhalation_rate=seated_act.exhalation_rate),
presence = interaction_intervals[0],
distance = 0.854,
),
mc.ShortRangeModel(
data_registry = data_registry,
expiration = short_range_expiration_distributions(data_registry)['Breathing'],
activity = models.Activity(inhalation_rate=heavy_exercise_act.inhalation_rate * 2,
exhalation_rate=heavy_exercise_act.inhalation_rate),
presence = interaction_intervals[1],
distance = 0.854,
),
)
e_model_2: models.ExposureModel = exposure_model_from_parameter(data_registry=data_registry, short_range_models=sr_models_activity, BR=2.5)
np.testing.assert_allclose(
2*e_model_1.deposited_exposure().mean(),
e_model_2.deposited_exposure().mean(), rtol=0.02
)