Skip to content

Commit

Permalink
Merge branch 'bugfix/short-range-f_inf' into 'master'
Browse files Browse the repository at this point in the history 
Fixed error with f_inf (short-range)

Closes #404

See merge request caimira/caimira!496
  • Loading branch information
@lrdossan
lrdossan committed Jul 15, 2024
2 parents e1a94b9 + bbb82d2 commit bad998c
Show file tree
Hide file tree
Showing 5 changed files with 220 additions and 73 deletions.
2 changes: 1 addition & 1 deletion caimira/apps/calculator/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@
# 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")

Expand Down
5 changes: 4 additions & 1 deletion caimira/apps/calculator/report_generator.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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')
Expand Down
140 changes: 80 additions & 60 deletions caimira/models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -687,13 +687,8 @@ def particle(self) -> Particle:

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

Expand Down Expand Up @@ -723,8 +718,9 @@ def particle(self) -> Particle:
@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.
Expand All @@ -733,14 +729,6 @@ def volume(d):
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):
Expand Down Expand Up @@ -889,9 +877,9 @@ def aerosols(self):

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")
Expand All @@ -900,12 +888,12 @@ def emission_rate_per_aerosol_per_person_when_present(self) -> _VectorisedFloat:
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.
"""
Expand Down Expand Up @@ -945,13 +933,13 @@ def aerosols(self):
@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):
Expand All @@ -974,17 +962,19 @@ def aerosols(self):
@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
Expand Down Expand Up @@ -1386,29 +1376,38 @@ def dilution_factor(self) -> _VectorisedFloat:
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:
Expand All @@ -1420,15 +1419,32 @@ def _normed_concentration(self, concentration_model: ConcentrationModel, time: f

# 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:
Expand Down Expand Up @@ -1462,37 +1478,41 @@ def extract_between_bounds(self, time1: float, time2: float) -> typing.Union[Non
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(
self, concentration_model: ConcentrationModel,
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(
Expand Down Expand Up @@ -1566,7 +1586,8 @@ def fun(x):
@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

Expand Down Expand Up @@ -1693,7 +1714,6 @@ def long_range_deposited_exposure_between_bounds(self, time1: float, time2: floa
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:
Expand All @@ -1709,8 +1729,7 @@ def deposited_exposure_between_bounds(self, time1: float, time2: float) -> _Vect
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()
Expand Down Expand Up @@ -1741,11 +1760,12 @@ def deposited_exposure_between_bounds(self, time1: float, time2: float) -> _Vect
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)

Expand Down
8 changes: 4 additions & 4 deletions caimira/tests/models/test_short_range_model.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down Expand Up @@ -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.],
]
)
Expand Down
Loading

0 comments on commit bad998c

Please sign in to comment.