Change all normalization factors and emission rates (when present) into 'per person' quantities (i.e. do not normalize anymore by the number of infected people) - to allow dynamic number of infected people. Changing the test accordingly, and the documentation. Adding a __post_init__ check in ExposureModel to forbid the use of dynamic number of exposed people, for now.
This commit is contained in:
Nicolas Mounet
Luis Aleixo
parent
e494151d78
commit
fe19cb5976
8 changed files with 98 additions and 113 deletions
|
|
@ -145,7 +145,7 @@ def calculate_report_data(form: FormData, model: models.ExposureModel) -> typing
|
|||
prob = np.array(model.infection_probability())
|
||||
prob_dist_count, prob_dist_bins = np.histogram(prob/100, bins=100, density=True)
|
||||
prob_probabilistic_exposure = np.array(model.total_probability_rule()).mean()
|
||||
er = np.array(model.concentration_model.infected.emission_rate_when_present()).mean()
|
||||
er = np.array(model.concentration_model.infected.emission_rate_per_person_when_present()).mean()
|
||||
exposed_occupants = model.exposed.number
|
||||
expected_new_cases = np.array(model.expected_new_cases()).mean()
|
||||
uncertainties_plot_src = img2base64(_figure2bytes(uncertainties_plot(model))) if form.conditional_probability_plot else None
|
||||
|
|
@ -166,12 +166,12 @@ def calculate_report_data(form: FormData, model: models.ExposureModel) -> typing
|
|||
"cumulative_doses": list(cumulative_doses),
|
||||
"long_range_cumulative_doses": list(long_range_cumulative_doses),
|
||||
"prob_inf": prob.mean(),
|
||||
"prob_inf_sd": np.std(prob),
|
||||
"prob_inf_sd": prob.std(),
|
||||
"prob_dist": list(prob),
|
||||
"prob_hist_count": list(prob_dist_count),
|
||||
"prob_hist_bins": list(prob_dist_bins),
|
||||
"prob_probabilistic_exposure": prob_probabilistic_exposure,
|
||||
"emission_rate": er,
|
||||
"emission_rate_per_person": er,
|
||||
"exposed_occupants": exposed_occupants,
|
||||
"expected_new_cases": expected_new_cases,
|
||||
"uncertainties_plot_src": uncertainties_plot_src,
|
||||
|
|
|
|||
|
|
@ -209,7 +209,7 @@ class ExposureModelResult(View):
|
|||
def update_textual_result(self, model: models.ExposureModel):
|
||||
lines = []
|
||||
P = np.array(model.infection_probability()).mean()
|
||||
lines.append(f'Emission rate (virus/hr): {np.round(model.concentration_model.infected.emission_rate_when_present(),0)}')
|
||||
lines.append(f'Emission rate per infected person (virus/hr): {np.round(model.concentration_model.infected.emission_rate_per_person_when_present(),0)}')
|
||||
lines.append(f'Probability of infection: {np.round(P, 0)}%')
|
||||
|
||||
lines.append(f'Number of exposed: {model.exposed.number}')
|
||||
|
|
|
|||
|
|
@ -26,7 +26,7 @@ provided the sample size is large enough. Example of the MC integration over the
|
|||
|
||||
It is important to distinguish between 1) Monte-Carlo random variables (which are vectorised independently on its diameter-dependence) and 2) numerical Monte-Carlo integration for the diameter-dependence.
|
||||
Since the integral of the diameter-dependent variables are solved when computing the dose -- :math:`\mathrm{vD^{total}}` -- while performing some of the intermediate calculations,
|
||||
we normalize the results by *dividing* by the Monte-Carlo variables that are diameter-independent, so that they are not considered in the Monte-Carlo integration (e.g. the **viral load** parameter, or the result of the :meth:`caimira.models.InfectedPopulation.emission_rate_per_aerosol_when_present` method).
|
||||
we normalize the results by *dividing* by the Monte-Carlo variables that are diameter-independent, so that they are not considered in the Monte-Carlo integration (e.g. the **viral load** parameter, or the result of the :meth:`caimira.models.InfectedPopulation.emission_rate_per_aerosol_per_person_when_present` method).
|
||||
|
||||
Expiration
|
||||
==========
|
||||
|
|
@ -62,9 +62,9 @@ Note that :math:`D_{\mathrm{max}}` value will differ, depending on the type of e
|
|||
|
||||
In the code, for a given Expiration, we use different methods to perform the calculations *step-by-step*:
|
||||
|
||||
1. Calculate the non aerosol-dependent quantities in the emission rate, which is the multiplication of the diameter-**independent** variables: :meth:`caimira.models.InfectedPopulation.emission_rate_per_aerosol_when_present`. This corresponds to the :math:`\mathrm{vl_{in}} \cdot \mathrm{BR_{k}}` part of the :math:`\mathrm{vR}(D)` equation.
|
||||
1. Calculate the non aerosol-dependent quantities in the emission rate per person infected, which is the multiplication of the diameter-**independent** variables: :meth:`caimira.models.InfectedPopulation.emission_rate_per_aerosol_per_person_when_present`. This corresponds to the :math:`\mathrm{vl_{in}} \cdot \mathrm{BR_{k}}` part of the :math:`\mathrm{vR}(D)` equation.
|
||||
2. Calculate the diameter-**dependent** variable :meth:`caimira.models.InfectedPopulation.aerosols`, which is the result of :math:`E_{c,j}(D) = N_p(D) \cdot V_p(D) \cdot (1 − η_\mathrm{out}(D))` (in mL/(m\ :sup:`3` \.µm)), with :math:`N_p(D)` being the product of the BLO distribution by the scaling factor :math:`cn`. Note that this result is not integrated over the diameters at this stage, thus the units are still *'per aerosol diameter'*.
|
||||
3. Calculate the full emission rate, which is the multiplication of the two previous methods, and corresponds to :math:`\mathrm{vR(D)}`: :meth:`caimira.models._PopulationWithVirus.emission_rate_when_present`.
|
||||
3. Calculate the full emission rate (per person infected), which is the multiplication of the two previous methods, and corresponds to :math:`\mathrm{vR(D)}`: :meth:`caimira.models._PopulationWithVirus.emission_rate_per_person_when_present`.
|
||||
|
||||
Note that the diameter-dependence is kept at this stage. Since other parameters downstream in code are also diameter-dependent, the Monte-Carlo integration over the aerosol sizes is computed at the level of the dose :math:`\mathrm{vD^{total}}`.
|
||||
In case one would like to have intermediate results for emission rate, perform the Monte-Carlo integration of :math:`E_{c, j}^{\mathrm{total}}` and compute :math:`\mathrm{vR^{total}} =\mathrm{vl_{in}} \cdot E_{c, j}^{\mathrm{total}} \cdot \mathrm{BR_k}`.
|
||||
|
|
@ -97,8 +97,8 @@ the code computes first a normalized version of the concentration, i.e. divided
|
|||
|
||||
To summarize, we can split the concentration in two different formulations:
|
||||
|
||||
* Normalized concentration :meth:`caimira.models._ConcentrationModelBase._normed_concentration`: :math:`\mathrm{C_\mathrm{LR, normed}}(t, D)` that computes the concentration without including the emission rate.
|
||||
* Concentration :meth:`caimira.models._ConcentrationModelBase.concentration` : :math:`C_{\mathrm{LR}}(t, D) = \mathrm{C_\mathrm{LR, normed}}(t, D) \cdot \mathrm{vR}(D)`, where :math:`\mathrm{vR}(D)` is the result of the :meth:`caimira.models._PopulationWithVirus.emission_rate_when_present` method.
|
||||
* Normalized concentration :meth:`caimira.models._ConcentrationModelBase._normed_concentration`: :math:`\mathrm{C_\mathrm{LR, normed}}(t, D)` that computes the concentration without including the emission rate per person infected.
|
||||
* Concentration :meth:`caimira.models._ConcentrationModelBase.concentration` : :math:`C_{\mathrm{LR}}(t, D) = \mathrm{C_\mathrm{LR, normed}}(t, D) \cdot \mathrm{vR}(D)`, where :math:`\mathrm{vR}(D)` is the result of the :meth:`caimira.models._PopulationWithVirus.emission_rate_per_person_when_present` method.
|
||||
|
||||
Note that in order to get the total concentration value in this stage, the final result should be averaged over the particle diameters (i.e. Monte-Carlo integration over diameters, see above).
|
||||
For the calculator app report, the total concentration (MC integral over the diameter) is performed only when generating the plot.
|
||||
|
|
@ -106,8 +106,8 @@ Otherwise, the diameter-dependence continues until we compute the inhaled dose i
|
|||
|
||||
The following methods calculate the integrated concentration between two times. They are mostly used when calculating the **dose**:
|
||||
|
||||
* :meth:`caimira.models._ConcentrationModelBase.normed_integrated_concentration`, :math:`\mathrm{C_\mathrm{normed}}(D)` that returns the integrated long-range concentration of viruses in the air, between any two times, normalized by the emission rate. Note that this method performs the integral between any two times of the previously mentioned :meth:`caimira.models._ConcentrationModelBase._normed_concentration` method.
|
||||
* :meth:`caimira.models._ConcentrationModelBase.integrated_concentration`, :math:`C(D)`, that returns the same result as the previous one, but multiplied by the emission rate.
|
||||
* :meth:`caimira.models._ConcentrationModelBase.normed_integrated_concentration`, :math:`\mathrm{C_\mathrm{normed}}(D)` that returns the integrated long-range concentration of viruses in the air, between any two times, normalized by the emission rate per person infected. Note that this method performs the integral between any two times of the previously mentioned :meth:`caimira.models._ConcentrationModelBase._normed_concentration` method.
|
||||
* :meth:`caimira.models._ConcentrationModelBase.integrated_concentration`, :math:`C(D)`, that returns the same result as the previous one, but multiplied by the emission rate (per person infected).
|
||||
|
||||
The integral over the exposure times is calculated directly in the class (integrated methods).
|
||||
|
||||
|
|
@ -221,14 +221,14 @@ Long-range approach
|
|||
*******************
|
||||
|
||||
Regarding the concentration part of the long-range exposure (concentration integrated over time, :math:`\int_{t1}^{t2}C_{\mathrm{LR}}(t, D)\;\mathrm{d}t`), the respective method is :meth:`caimira.models.ExposureModel._long_range_normed_exposure_between_bounds`,
|
||||
which uses the long-range exposure (concentration) between two bounds (time1 and time2), normalized by the emission rate of the infected population, calculated from :meth:`caimira.models._ConcentrationModelBase.normed_integrated_concentration`.
|
||||
which uses the long-range exposure (concentration) between two bounds (time1 and time2), normalized by the emission rate of the infected population (per person infected), calculated from :meth:`caimira.models._ConcentrationModelBase.normed_integrated_concentration`.
|
||||
The former method filters out the given bounds considering the breaks through the day (i.e. the time intervals during which there is no exposition to the virus) and retrieves the integrated long-range concentration of viruses in the air between any two times.
|
||||
|
||||
After the calculations of the integrated concentration over the time, in order to calculate the final dose, we have to compute the remaining factors in the above equation.
|
||||
Note that the **Monte-Carlo integration over the diameters is performed at this stage**, where all the diameter-dependent parameters are grouped together to calculate the final average (:code:`np.mean()`).
|
||||
|
||||
Since, in the previous chapters, the quantities where normalised by the emission rate, one will need to re-incorporate it in the equations before performing the MC integrations over :math:`D`.
|
||||
For that we need to split :math:`\mathrm{vR}(D)` (:meth:`caimira.models._PopulationWithVirus.emission_rate_when_present`) in diameter-dependent and diameter-independent quantities:
|
||||
Since, in the previous chapters, the quantities where normalised by the emission rate per person infected, one will need to re-incorporate it in the equations before performing the MC integrations over :math:`D`.
|
||||
For that we need to split :math:`\mathrm{vR}(D)` (:meth:`caimira.models._PopulationWithVirus.emission_rate_per_person_when_present`) in diameter-dependent and diameter-independent quantities:
|
||||
|
||||
:math:`\mathrm{vR}(D) = \mathrm{vR}(D-\mathrm{dependent}) \times \mathrm{vR}(D-\mathrm{independent})`
|
||||
|
||||
|
|
@ -236,7 +236,7 @@ with
|
|||
|
||||
:math:`\mathrm{vR}(D-\mathrm{dependent}) = \mathrm{cn} \cdot V_p(D) \cdot (1 − \mathrm{η_{out}}(D))` - :meth:`caimira.models.InfectedPopulation.aerosols`
|
||||
|
||||
:math:`\mathrm{vR}(D-\mathrm{independent}) = \mathrm{vl_{in}} \cdot \mathrm{BR_{k}}` - :meth:`caimira.models.InfectedPopulation.emission_rate_per_aerosol_when_present`
|
||||
:math:`\mathrm{vR}(D-\mathrm{independent}) = \mathrm{vl_{in}} \cdot \mathrm{BR_{k}}` - :meth:`caimira.models.InfectedPopulation.emission_rate_per_aerosol_per_person_when_present`
|
||||
|
||||
|
||||
In other words, in the code the procedure is the following (all performed in :meth:`caimira.models.ExposureModel.long_range_deposited_exposure_between_bounds` method):
|
||||
|
|
@ -248,7 +248,7 @@ In other words, in the code the procedure is the following (all performed in :me
|
|||
* in order to complete the equation, multiply by the remaining diameter-independent variables in :math:`\mathrm{vD}` to obtain the total value: :math:`\mathrm{vD^{total}} = \mathrm{vD_{emission-rate}} \cdot \mathrm{BR}_{\mathrm{k}} \cdot (1-\eta_{\mathrm{in}}) \cdot f_{\mathrm{inf}}`;
|
||||
* in the end, the dose is a vectorized float used in the probability of infection formula.
|
||||
|
||||
**Note**: The aerosol volume concentration (*aerosols*) is introduced because the integrated concentration over the time was previously normalized by the emission rate.
|
||||
**Note**: The aerosol volume concentration (*aerosols*) is introduced because the integrated concentration over the time was previously normalized by the emission rate (per person).
|
||||
Here, to calculate the integral over the diameters we also need to consider the diameter-dependent variables that are on the emission rate, represented by the aerosol volume concentration which depends on the diameter and on the mask type:
|
||||
|
||||
:math:`\mathrm{aerosols} = \mathrm{cn} \cdot V_p(D) \cdot (1 − \mathrm{η_{out}}(D))` .
|
||||
|
|
@ -339,4 +339,4 @@ REFERENCES
|
|||
==========
|
||||
|
||||
.. [1] Jia, Wei, et al. "Exposure and respiratory infection risk via the short-range airborne route." Building and environment 219 (2022): 109166. `doi.org/10.1016/j.buildenv.2022.109166 <https://doi.org/10.1016/j.buildenv.2022.109166>`_
|
||||
.. [2] Henriques, Andre, et al. "Modelling airborne transmission of SARS-CoV-2 using CARA: risk assessment for enclosed spaces." Interface Focus 12.2 (2022): 20210076. `doi.org/10.1098/rsfs.2021.0076 <https://doi.org/10.1098/rsfs.2021.0076>`_
|
||||
.. [2] Henriques, Andre, et al. "Modelling airborne transmission of SARS-CoV-2 using CARA: risk assessment for enclosed spaces." Interface Focus 12.2 (2022): 20210076. `doi.org/10.1098/rsfs.2021.0076 <https://doi.org/10.1098/rsfs.2021.0076>`_
|
||||
|
|
|
|||
|
|
@ -819,15 +819,6 @@ class Population:
|
|||
else:
|
||||
if self.presence is not None:
|
||||
raise TypeError(f'The presence argument must be None for a IntPiecewiseConstant number.')
|
||||
|
||||
def _last_presence_state(self, time: float):
|
||||
"""
|
||||
Find the most recent/previous presence state change.
|
||||
"""
|
||||
times = set(self.number.transition_times) # type: ignore
|
||||
t_index: int = np.searchsorted(sorted(times), time) # type: ignore
|
||||
t_index = max([t_index - 1, 0])
|
||||
return sorted(times)[t_index]
|
||||
|
||||
def person_present(self, time: float):
|
||||
# Allow back-compatibility
|
||||
|
|
@ -835,25 +826,15 @@ class Population:
|
|||
return self.presence.triggered(time)
|
||||
elif isinstance(self.number, IntPiecewiseContant):
|
||||
return self.number.value(time) != 0
|
||||
|
||||
def people_present(self, time: typing.Optional[float] = None):
|
||||
|
||||
def people_present(self, time: float):
|
||||
# Allow back-compatibility
|
||||
if isinstance(self.number, int) and isinstance(self.presence, Interval):
|
||||
return self.number
|
||||
|
||||
elif isinstance(self.number, IntPiecewiseContant):
|
||||
# For a given scenario, time falls within a break,
|
||||
# we should get the previous state occupancy profile.
|
||||
if time and not self.person_present(time):
|
||||
if time <= self.number.transition_times[0]:
|
||||
number = self.number.values[0]
|
||||
elif time >= self.number.transition_times[-1]:
|
||||
number = self.number.values[-1]
|
||||
else:
|
||||
number = int(self.number.value(self._last_presence_state(time)))
|
||||
else:
|
||||
number = int(self.number.value(time))
|
||||
return number
|
||||
if isinstance(self.number, int):
|
||||
return self.number * self.person_present(time)
|
||||
|
||||
else:
|
||||
return int(self.number.value(time))
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class _PopulationWithVirus(Population):
|
||||
|
|
@ -874,22 +855,22 @@ class _PopulationWithVirus(Population):
|
|||
"""
|
||||
raise NotImplementedError("Subclass must implement")
|
||||
|
||||
def emission_rate_per_aerosol_when_present(self, time: typing.Optional[float] = None) -> _VectorisedFloat:
|
||||
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).
|
||||
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.
|
||||
It should not be a function of time.
|
||||
"""
|
||||
raise NotImplementedError("Subclass must implement")
|
||||
|
||||
@method_cache
|
||||
def emission_rate_when_present(self, time: typing.Optional[float] = None) -> _VectorisedFloat:
|
||||
def emission_rate_per_person_when_present(self) -> _VectorisedFloat:
|
||||
"""
|
||||
The emission rate if the infected population is present
|
||||
The emission rate if the infected population is present, per person
|
||||
(in virions / h).
|
||||
"""
|
||||
return (self.emission_rate_per_aerosol_when_present(time) *
|
||||
return (self.emission_rate_per_aerosol_per_person_when_present() *
|
||||
self.aerosols())
|
||||
|
||||
def emission_rate(self, time) -> _VectorisedFloat:
|
||||
|
|
@ -906,7 +887,7 @@ class _PopulationWithVirus(Population):
|
|||
# itself a function of time. Any change in rate must be accompanied
|
||||
# with a declaration of state change time, as is the case for things
|
||||
# like Ventilation.
|
||||
return self.emission_rate_when_present(time)
|
||||
return self.emission_rate_per_person_when_present() * self.people_present(time)
|
||||
|
||||
@property
|
||||
def particle(self) -> Particle:
|
||||
|
|
@ -930,14 +911,14 @@ class EmittingPopulation(_PopulationWithVirus):
|
|||
return 1.
|
||||
|
||||
@method_cache
|
||||
def emission_rate_per_aerosol_when_present(self, time: typing.Optional[float] = None) -> _VectorisedFloat:
|
||||
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).
|
||||
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.
|
||||
It should not be a function of time.
|
||||
"""
|
||||
return self.known_individual_emission_rate * self.people_present(time)
|
||||
return self.known_individual_emission_rate
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
|
|
@ -959,7 +940,7 @@ class InfectedPopulation(_PopulationWithVirus):
|
|||
return self.expiration.aerosols(self.mask)
|
||||
|
||||
@method_cache
|
||||
def emission_rate_per_aerosol_when_present(self, time: typing.Optional[float] = None) -> _VectorisedFloat:
|
||||
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).
|
||||
|
|
@ -972,7 +953,7 @@ class InfectedPopulation(_PopulationWithVirus):
|
|||
ER = (self.virus.viral_load_in_sputum *
|
||||
self.activity.exhalation_rate *
|
||||
10 ** 6)
|
||||
return ER * self.people_present(time)
|
||||
return ER
|
||||
|
||||
@property
|
||||
def particle(self) -> Particle:
|
||||
|
|
@ -1040,7 +1021,7 @@ class _ConcentrationModelBase:
|
|||
"""
|
||||
return 0.
|
||||
|
||||
def normalization_factor(self, time: typing.Optional[float] = None) -> _VectorisedFloat:
|
||||
def normalization_factor(self) -> _VectorisedFloat:
|
||||
"""
|
||||
Normalization factor (in the same unit as the concentration).
|
||||
This factor is applied to the normalized concentration only
|
||||
|
|
@ -1059,13 +1040,14 @@ class _ConcentrationModelBase:
|
|||
can be put back in front of the concentration after the time
|
||||
dependence has been solved for.
|
||||
"""
|
||||
if not self.population.person_present(time):
|
||||
return self.min_background_concentration()/self.normalization_factor(time)
|
||||
#if not self.population.person_present(time):
|
||||
# return self.min_background_concentration()/self.normalization_factor()
|
||||
|
||||
V = self.room.volume
|
||||
RR = self.removal_rate(time)
|
||||
|
||||
return (1. / (RR * V) + self.min_background_concentration()/
|
||||
self.normalization_factor(time))
|
||||
return (self.population.people_present(time) / (RR * V) +
|
||||
self.min_background_concentration()/self.normalization_factor())
|
||||
|
||||
@method_cache
|
||||
def state_change_times(self) -> typing.List[float]:
|
||||
|
|
@ -1143,7 +1125,7 @@ class _ConcentrationModelBase:
|
|||
# The model always starts at t=0, but we avoid running concentration calculations
|
||||
# before the first presence as an optimisation.
|
||||
if time <= self._first_presence_time():
|
||||
return self.min_background_concentration()/self.normalization_factor(time)
|
||||
return self.min_background_concentration()/self.normalization_factor()
|
||||
|
||||
next_state_change_time = self._next_state_change(time)
|
||||
|
||||
|
|
@ -1156,12 +1138,12 @@ class _ConcentrationModelBase:
|
|||
conc_limit = 0.
|
||||
|
||||
t_last_state_change = self.last_state_change(time)
|
||||
conc_at_last_state_change = self._normed_concentration_cached(t_last_state_change)*self.normalization_factor(t_last_state_change)
|
||||
conc_at_last_state_change = self._normed_concentration_cached(t_last_state_change)
|
||||
|
||||
delta_time = time - t_last_state_change
|
||||
fac = np.exp(-RR * delta_time)
|
||||
|
||||
return conc_limit * (1 - fac) + conc_at_last_state_change/self.normalization_factor(time) * fac
|
||||
return conc_limit * (1 - fac) + conc_at_last_state_change * fac
|
||||
|
||||
def concentration(self, time: float) -> _VectorisedFloat:
|
||||
"""
|
||||
|
|
@ -1172,7 +1154,7 @@ class _ConcentrationModelBase:
|
|||
to this method.
|
||||
"""
|
||||
return (self._normed_concentration_cached(time) *
|
||||
self.normalization_factor(time))
|
||||
self.normalization_factor())
|
||||
|
||||
@method_cache
|
||||
def normed_integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat:
|
||||
|
|
@ -1209,7 +1191,7 @@ class _ConcentrationModelBase:
|
|||
Get the integrated concentration of viruses in the air between the times start and stop.
|
||||
"""
|
||||
return (self.normed_integrated_concentration(start, stop) *
|
||||
self.normalization_factor(start))
|
||||
self.normalization_factor())
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
|
|
@ -1234,9 +1216,9 @@ class ConcentrationModel(_ConcentrationModelBase):
|
|||
def virus(self) -> Virus:
|
||||
return self.infected.virus
|
||||
|
||||
def normalization_factor(self, time: typing.Optional[float] = None) -> _VectorisedFloat:
|
||||
def normalization_factor(self) -> _VectorisedFloat:
|
||||
# we normalize by the emission rate
|
||||
return self.infected.emission_rate_when_present(time)
|
||||
return self.infected.emission_rate_per_person_when_present()
|
||||
|
||||
def removal_rate(self, time: float) -> _VectorisedFloat:
|
||||
# Equilibrium velocity of particle motion toward the floor
|
||||
|
|
@ -1284,11 +1266,10 @@ class CO2ConcentrationModel(_ConcentrationModelBase):
|
|||
"""
|
||||
return self.CO2_atmosphere_concentration
|
||||
|
||||
def normalization_factor(self, time: typing.Optional[float] = None) -> _VectorisedFloat:
|
||||
# normalization by the CO2 exhaled.
|
||||
def normalization_factor(self) -> _VectorisedFloat:
|
||||
# normalization by the CO2 exhaled per person.
|
||||
# CO2 concentration given in ppm, hence the 1e6 factor.
|
||||
return (1e6*self.population.people_present(time)
|
||||
*self.population.activity.exhalation_rate
|
||||
return (1e6*self.population.activity.exhalation_rate
|
||||
*self.CO2_fraction_exhaled)
|
||||
|
||||
|
||||
|
|
@ -1515,8 +1496,11 @@ class ExposureModel:
|
|||
and not (
|
||||
all(np.isscalar(c_model.virus.decay_constant(c_model.room.humidity, c_model.room.inside_temp.value(time)) +
|
||||
c_model.ventilation.air_exchange(c_model.room, time)) for time in c_model.state_change_times()))):
|
||||
raise ValueError("If the diameter is an array, none of the ventilation parameters "
|
||||
"or virus decay constant can be arrays at the same time.")
|
||||
raise ValueError("If the diameter is an array, none of the ventilation parameters "
|
||||
"or virus decay constant can be arrays at the same time.")
|
||||
if not isinstance(self.exposed.number, int):
|
||||
raise NotImplementedError("Cannot use dynamic occupancy for"
|
||||
" the exposed population")
|
||||
|
||||
def long_range_fraction_deposited(self) -> _VectorisedFloat:
|
||||
"""
|
||||
|
|
@ -1533,21 +1517,19 @@ class ExposureModel:
|
|||
by the emission rate of the infected population
|
||||
"""
|
||||
exposure = 0.
|
||||
if isinstance(self.exposed.number, int) and isinstance(self.exposed.presence, Interval):
|
||||
# TODO: Implement the case for the IntPiecewiseConstant
|
||||
for start, stop in self.exposed.presence.boundaries():
|
||||
if stop < time1:
|
||||
continue
|
||||
elif start > time2:
|
||||
break
|
||||
elif start <= time1 and time2<= stop:
|
||||
exposure += self.concentration_model.normed_integrated_concentration(time1, time2)
|
||||
elif start <= time1 and stop < time2:
|
||||
exposure += self.concentration_model.normed_integrated_concentration(time1, stop)
|
||||
elif time1 < start and time2 <= stop:
|
||||
exposure += self.concentration_model.normed_integrated_concentration(start, time2)
|
||||
elif time1 <= start and stop < time2:
|
||||
exposure += self.concentration_model.normed_integrated_concentration(start, stop)
|
||||
for start, stop in self.exposed.presence.boundaries():
|
||||
if stop < time1:
|
||||
continue
|
||||
elif start > time2:
|
||||
break
|
||||
elif start <= time1 and time2<= stop:
|
||||
exposure += self.concentration_model.normed_integrated_concentration(time1, time2)
|
||||
elif start <= time1 and stop < time2:
|
||||
exposure += self.concentration_model.normed_integrated_concentration(time1, stop)
|
||||
elif time1 < start and time2 <= stop:
|
||||
exposure += self.concentration_model.normed_integrated_concentration(start, time2)
|
||||
elif time1 <= start and stop < time2:
|
||||
exposure += self.concentration_model.normed_integrated_concentration(start, stop)
|
||||
return exposure
|
||||
|
||||
def concentration(self, time: float) -> _VectorisedFloat:
|
||||
|
|
@ -1565,7 +1547,8 @@ class ExposureModel:
|
|||
def long_range_deposited_exposure_between_bounds(self, time1: float, time2: float) -> _VectorisedFloat:
|
||||
deposited_exposure = 0.
|
||||
|
||||
emission_rate_per_aerosol = self.concentration_model.infected.emission_rate_per_aerosol_when_present()
|
||||
emission_rate_per_aerosol_per_person = \
|
||||
self.concentration_model.infected.emission_rate_per_aerosol_per_person_when_present()
|
||||
aerosols = self.concentration_model.infected.aerosols()
|
||||
f_inf = self.concentration_model.infected.fraction_of_infectious_virus()
|
||||
fdep = self.long_range_fraction_deposited()
|
||||
|
|
@ -1585,10 +1568,11 @@ class ExposureModel:
|
|||
# one should not take any mean at this stage.
|
||||
dep_exposure_integrated = self._long_range_normed_exposure_between_bounds(time1, time2)*aerosols*fdep
|
||||
|
||||
# Then we multiply by the diameter-independent quantity emission_rate_per_aerosol,
|
||||
# Then we multiply by the diameter-independent quantity emission_rate_per_aerosol_per_person,
|
||||
# and parameters of the vD equation (i.e. BR_k and n_in).
|
||||
deposited_exposure += (dep_exposure_integrated * emission_rate_per_aerosol *
|
||||
self.exposed.activity.inhalation_rate *
|
||||
deposited_exposure += (dep_exposure_integrated *
|
||||
emission_rate_per_aerosol_per_person *
|
||||
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.
|
||||
|
|
@ -1655,10 +1639,8 @@ class ExposureModel:
|
|||
The number of virus per m^3 deposited on the respiratory tract.
|
||||
"""
|
||||
deposited_exposure: _VectorisedFloat = 0.0
|
||||
if isinstance(self.exposed.number, int) and isinstance(self.exposed.presence, Interval):
|
||||
# TODO: Implement the case for the IntPiecewiseConstant
|
||||
for start, stop in self.exposed.presence.boundaries():
|
||||
deposited_exposure += self.deposited_exposure_between_bounds(start, stop)
|
||||
for start, stop in self.exposed.presence.boundaries():
|
||||
deposited_exposure += self.deposited_exposure_between_bounds(start, stop)
|
||||
|
||||
return deposited_exposure * self.repeats
|
||||
|
||||
|
|
@ -1676,11 +1658,13 @@ class ExposureModel:
|
|||
def total_probability_rule(self) -> _VectorisedFloat:
|
||||
if (self.geographical_data.geographic_population != 0 and self.geographical_data.geographic_cases != 0):
|
||||
sum_probability = 0.0
|
||||
if not isinstance(self.concentration_model.infected.number, int):
|
||||
raise NotImplementedError("Cannot compute total probability "
|
||||
"(including incidence rate) with dynamic occupancy")
|
||||
|
||||
# Create an equivalent exposure model but changing the number of infected cases.
|
||||
if isinstance(self.concentration_model.infected.number, int) and isinstance(self.exposed.number, int):
|
||||
# TODO: Implement the case for the IntPiecewiseConstant
|
||||
total_people = self.concentration_model.infected.number + self.exposed.number
|
||||
max_num_infected = (total_people if total_people < 10 else 10)
|
||||
total_people = self.concentration_model.infected.number + self.exposed.number
|
||||
max_num_infected = (total_people if total_people < 10 else 10)
|
||||
# The influence of a higher number of simultainious infected people (> 4 - 5) yields an almost negligible contirbution to the total probability.
|
||||
# To be on the safe side, a hard coded limit with a safety margin of 2x was set.
|
||||
# Therefore we decided a hard limit of 10 infected people.
|
||||
|
|
|
|||
|
|
@ -33,7 +33,7 @@ class KnownConcentrationModelBase(models._ConcentrationModelBase):
|
|||
def min_background_concentration(self) -> float:
|
||||
return self.known_min_background_concentration
|
||||
|
||||
def normalization_factor(self, time: typing.Optional[float] = None) -> float:
|
||||
def normalization_factor(self) -> float:
|
||||
return self.known_normalization_factor
|
||||
|
||||
|
||||
|
|
@ -107,7 +107,7 @@ def simple_conc_model():
|
|||
@pytest.fixture
|
||||
def dummy_population(simple_conc_model) -> models.Population:
|
||||
return models.Population(
|
||||
number=10,
|
||||
number=1,
|
||||
presence=simple_conc_model.infected.presence,
|
||||
mask=models.Mask.types['Type I'],
|
||||
activity=models.Activity.types['Seated'],
|
||||
|
|
@ -268,7 +268,7 @@ def test_zero_ventilation_rate(
|
|||
ventilation = simple_conc_model.ventilation,
|
||||
known_population = dummy_population,
|
||||
known_removal_rate = known_removal_rate,
|
||||
known_normalization_factor=10.,
|
||||
known_normalization_factor=1.,
|
||||
known_min_background_concentration = known_min_background_concentration)
|
||||
|
||||
normed_concentration = known_conc_model.concentration(1)
|
||||
|
|
|
|||
|
|
@ -175,4 +175,4 @@ def test_dynamic_dose(full_exposure_model):
|
|||
dynamic_exposure = dynamic_infected.deposited_exposure()
|
||||
static_exposure = np.sum([model.deposited_exposure() for model in (single_infected, two_infected, three_infected)])
|
||||
|
||||
npt.assert_almost_equal(dynamic_exposure, static_exposure)
|
||||
npt.assert_almost_equal(dynamic_exposure, static_exposure)
|
||||
|
|
|
|||
|
|
@ -24,7 +24,7 @@ class KnownNormedconcentration(models.ConcentrationModel):
|
|||
return 1.e50
|
||||
|
||||
def _normed_concentration_limit(self, time: float) -> models._VectorisedFloat:
|
||||
return self.normed_concentration_function(time)
|
||||
return self.normed_concentration_function(time) * self.infected.number
|
||||
|
||||
def state_change_times(self):
|
||||
return [0., 24.]
|
||||
|
|
@ -33,7 +33,7 @@ class KnownNormedconcentration(models.ConcentrationModel):
|
|||
return 24.
|
||||
|
||||
def _normed_concentration(self, time: float) -> models._VectorisedFloat: # noqa
|
||||
return self.normed_concentration_function(time)
|
||||
return self.normed_concentration_function(time) * self.infected.number
|
||||
|
||||
|
||||
halftime = models.PeriodicInterval(120, 60)
|
||||
|
|
@ -67,7 +67,8 @@ def known_concentrations(func):
|
|||
expiration=models.Expiration.types['Speaking'],
|
||||
host_immunity=0.,
|
||||
)
|
||||
normed_func = lambda x: func(x) / dummy_infected_population.emission_rate_when_present()
|
||||
normed_func = lambda x: (func(x) /
|
||||
dummy_infected_population.emission_rate_per_person_when_present())
|
||||
return KnownNormedconcentration(dummy_room, dummy_ventilation,
|
||||
dummy_infected_population, 0.3, normed_func)
|
||||
|
||||
|
|
|
|||
|
|
@ -313,19 +313,19 @@ def waiting_room_mc():
|
|||
|
||||
@retry(tries=10)
|
||||
@pytest.mark.parametrize(
|
||||
"mc_model, expected_pi, expected_new_cases, expected_dose, expected_ER",
|
||||
"mc_model, expected_pi, expected_new_cases, expected_dose, expected_ER_per_person",
|
||||
[
|
||||
["shared_office_mc", 5.38, 0.16, 3.350, 1056],
|
||||
["classroom_mc", 8.21, 1.56, 11.356, 7416],
|
||||
["ski_cabin_mc", 12.92, 0.39, 21.796, 10231],
|
||||
["skagit_chorale_mc",61.01, 36.53, 84.730, 190422],
|
||||
["bus_ride_mc", 10.59, 7.06, 6.650, 5419],
|
||||
["gym_mc", 0.52, 0.14, 0.249, 1450],
|
||||
["gym_mc", 0.52, 0.14, 0.249, 1450/2.], # there are two infected in this case
|
||||
["waiting_room_mc", 1.53, 0.21, 0.844, 929],
|
||||
]
|
||||
)
|
||||
def test_report_models(mc_model, expected_pi, expected_new_cases,
|
||||
expected_dose, expected_ER, request):
|
||||
expected_dose, expected_ER_per_person, request):
|
||||
mc_model = request.getfixturevalue(mc_model)
|
||||
exposure_model = mc_model.build_model(size=SAMPLE_SIZE)
|
||||
npt.assert_allclose(exposure_model.infection_probability().mean(),
|
||||
|
|
@ -335,8 +335,8 @@ def test_report_models(mc_model, expected_pi, expected_new_cases,
|
|||
npt.assert_allclose(exposure_model.deposited_exposure().mean(),
|
||||
expected_dose, rtol=TOLERANCE)
|
||||
npt.assert_allclose(
|
||||
exposure_model.concentration_model.infected.emission_rate_when_present().mean(),
|
||||
expected_ER, rtol=TOLERANCE)
|
||||
exposure_model.concentration_model.infected.emission_rate_per_person_when_present().mean(),
|
||||
expected_ER_per_person, rtol=TOLERANCE)
|
||||
|
||||
|
||||
@retry(tries=10)
|
||||
|
|
@ -397,5 +397,5 @@ def test_small_shared_office_Geneva(mask_type, month, expected_pi,
|
|||
npt.assert_allclose(exposure_model.deposited_exposure().mean(),
|
||||
expected_dose, rtol=TOLERANCE)
|
||||
npt.assert_allclose(
|
||||
exposure_model.concentration_model.infected.emission_rate_when_present().mean(),
|
||||
exposure_model.concentration_model.infected.emission_rate_per_person_when_present().mean(),
|
||||
expected_ER, rtol=TOLERANCE)
|
||||
|
|
|
|||
Loading…
Reference in a new issue