diff --git a/caimira/apps/calculator/__init__.py b/caimira/apps/calculator/__init__.py index 50790196..2f42bd0d 100644 --- a/caimira/apps/calculator/__init__.py +++ b/caimira/apps/calculator/__init__.py @@ -35,7 +35,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.5" +__version__ = "4.6" class BaseRequestHandler(RequestHandler): diff --git a/caimira/apps/expert.py b/caimira/apps/expert.py index 6d489697..f48d7aaa 100644 --- a/caimira/apps/expert.py +++ b/caimira/apps/expert.py @@ -1091,4 +1091,4 @@ def models_start_end(models: typing.Sequence[models.ExposureModel]) -> typing.Tu """ infected_start = min(model.concentration_model.infected.presence.boundaries()[0][0] for model in models) infected_finish = min(model.concentration_model.infected.presence.boundaries()[-1][1] for model in models) - return infected_start, infected_finish \ No newline at end of file + return infected_start, infected_finish diff --git a/caimira/docs/UML-CAiMIRA.png b/caimira/docs/UML-CAiMIRA.png index 35fd4fb3..6b91e30b 100644 Binary files a/caimira/docs/UML-CAiMIRA.png and b/caimira/docs/UML-CAiMIRA.png differ diff --git a/caimira/docs/full_diameter_dependence.rst b/caimira/docs/full_diameter_dependence.rst index 3fbea5ba..1e684ebb 100644 --- a/caimira/docs/full_diameter_dependence.rst +++ b/caimira/docs/full_diameter_dependence.rst @@ -69,14 +69,18 @@ In the code, for a given Expiration, we use different methods to perform the cal 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}`. -Concentration - C(t, D) -======================= +Virus Concentration - C(t, D) +============================= The estimate of the concentration of virus-laden particles in a given room is based on a two-box exposure model: * **Box 1** - long-range exposure: also known as the *background* concentration, corresponds to the exposure of airborne virions where the susceptible (exposed) host is more than 2 m away from the infected host(s), considering the result of a mass balance equation between the emission rate of the infected host(s) and the removal rates from the environmental/virological characteristics. * **Box 2** - short-range exposure: also known as the *exhaled jet* concentration in close-proximity, corresponds to the exposure of airborne virions where the susceptible (exposed) host is distanced between 0.5 and 2 m from an infected host, considering the result of a two-stage exhaled jet model. +Note that most of the methods used to calculate the concentration are defined in the superclass :meth:`caimira.models._ConcentrationModelBase`, while the specific methods for the long-range virus concentration are part of the subclass :meth:`caimira.models.ConcentrationModel`. +The specific removal rate, minimum background concentration and normalization factors will depend on what concentration is being calculated (e.g. viral concentration or CO\ :sub:`2` concentration) and are respectively defined in :meth:`caimira.models._ConcentrationModelBase.removal_rate`, +:meth:`caimira.models._ConcentrationModelBase.min_background_concentration` and :meth:`caimira.models._ConcentrationModelBase.normalization_factor`. + Long-range approach ******************* @@ -84,8 +88,8 @@ The long-range concentration of virus-laden aerosols of a given size :math:`D`, :math:`C_{\mathrm{LR}}(t, D)=\frac{\mathrm{vR}(D) \cdot N_{\mathrm{inf}}}{\lambda_{\mathrm{vRR}}(D) \cdot V_r}-\left (\frac{\mathrm{vR}(D) \cdot N_{\mathrm{inf}}}{\lambda_{\mathrm{vRR}}(D) \cdot V_r}-C_0(D) \right )e^{-\lambda_{\mathrm{vRR}}(D)t}` , -and computed, as a function of the exposure time and particle diameter, in the :meth:`caimira.models.ConcentrationModel.concentration` method. -The long-range concentration, integrated over the exposure time (in piecewise constant steps), :math:`C(D)`, is given by :meth:`caimira.models.ConcentrationModel.integrated_concentration`. +and computed, as a function of the exposure time and particle diameter, in the :meth:`caimira.models._ConcentrationModelBase.concentration` method. +The long-range concentration, integrated over the exposure time (in piecewise constant steps), :math:`C(D)`, is given by :meth:`caimira.models._ConcentrationModelBase.integrated_concentration`. In the :math:`C_{\mathrm{LR}}(t, D)` equation above, the **emission rate** - :math:`\mathrm{vR}(D)` - and the **viral removal rate** - :math:`\lambda_{\mathrm{vRR}}(D)`, :meth:`caimira.models.ConcentrationModel.infectious_virus_removal_rate` - are both diameter-dependent. One can show that the resulting concentration is always proportional to the emission rate :math:`\mathrm{vR}(D)`. Hence, for computational speed-up purposes @@ -93,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.ConcentrationModel._normed_concentration`: :math:`\mathrm{C_\mathrm{LR, normed}}(t, D)` that computes the concentration without including the emission rate. -* Concentration :meth:`caimira.models.ConcentrationModel.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. +* 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. 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. @@ -102,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.ConcentrationModel.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.ConcentrationModel._normed_concentration` method. -* :meth:`caimira.models.ConcentrationModel.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. 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. The integral over the exposure times is calculated directly in the class (integrated methods). @@ -115,7 +119,7 @@ The short-range concentration is the result of a two-stage exhaled jet model dev :math:`C_{\mathrm{SR}}(t, D) = C_{\mathrm{LR}} (t, D) + \frac{1}{S({x})} \cdot (C_{0, \mathrm{SR}}(D) - C_{\mathrm{LR}, 100μm}(t, D))` , where :math:`S(x)` is the dilution factor due to jet dynamics, as a function of the interpersonal distance :math:`x` and :math:`C_{0, \mathrm{SR}}(D)` corresponds to the initial concentration of virions at the mouth/nose outlet during exhalation. -:math:`C_{\mathrm{LR}, 100μm}(t, D)` is the long-range concentration, calculated in :meth:`caimira.models.ConcentrationModel.concentration` method but **interpolated** to the diameter range used for close-proximity (from 0 to 100μm). +:math:`C_{\mathrm{LR}, 100μm}(t, D)` is the long-range concentration, calculated in :meth:`caimira.models._ConcentrationModelBase.concentration` method but **interpolated** to the diameter range used for close-proximity (from 0 to 100μm). Note that :math:`C_{0, \mathrm{SR}}(D)` is constant over time, hence only dependent on the particle diameter distribution. For code simplification, we split the :math:`C_{\mathrm{SR}}(t, D)` equation into two components: @@ -182,9 +186,9 @@ To summarize, in the code, :math:`C_{\mathrm{SR}}(t, D)` is computed as follows: * calculate the `dilution_factor` - :math:`S({x})` - in the method :meth:`caimira.models.ShortRangeModel.dilution_factor`, with the distance :math:`x` as a random variable (log normal distribution in :meth:`caimira.monte_carlo.data.short_range_distances`) * compute :math:`\frac{1}{S({x})} \cdot (C_{0, \mathrm{SR}}(D) - C_{\mathrm{LR}, 100\mathrm{μm}}(t, D))` in method :meth:`caimira.models.ShortRangeModel.normed_concentration`, * multiply by the diameter-independent parameter, viral load, in method :meth:`caimira.models.ShortRangeModel.short_range_concentration` -* complete the equation of :math:`C_{\mathrm{SR}}(t, D)` by adding the long-range concentration from the :meth:`caimira.models.ConcentrationModel.concentration` (all integrated over :math:`D`), returning the final short-range concentration value for a given time and expiration activity. This is done at the level of the Exposure Model (:meth:`caimira.models.ExposureModel.concentration`). +* complete the equation of :math:`C_{\mathrm{SR}}(t, D)` by adding the long-range concentration from the :meth:`caimira.models._ConcentrationModelBase.concentration` (all integrated over :math:`D`), returning the final short-range concentration value for a given time and expiration activity. This is done at the level of the Exposure Model (:meth:`caimira.models.ExposureModel.concentration`). -Note that :meth:`caimira.models.ShortRangeModel._normed_concentration` method is different from :meth:`caimira.models.ConcentrationModel._normed_concentration` and :meth:`caimira.models.ConcentrationModel.concentration` differs from :meth:`caimira.models.ExposureModel.concentration`. +Note that :meth:`caimira.models.ShortRangeModel._normed_concentration` method is different from :meth:`caimira.models._ConcentrationModelBase._normed_concentration` and :meth:`caimira.models._ConcentrationModelBase.concentration` differs from :meth:`caimira.models.ExposureModel.concentration`. Unless one is computing the mean concentration values (e.g. for the plots in the report), the diameter-dependence is kept at this stage. Since other parameters downstream in the code are also diameter-dependent, the Monte-Carlo integration over the particle sizes is computed at the level of the dose :math:`\mathrm{vD^{total}}`. In case one would like to have intermediate results for the initial short-range concentration, this is done at the :class:`caimira.models.ExposureModel` class level. @@ -217,7 +221,7 @@ 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.ConcentrationModel.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, 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. @@ -304,6 +308,20 @@ If short-range interactions exist: the long-range component is added to the alre If the are no short-range interactions: the short-range component (`deposited_exposure`) is zero, hence the result is equal solely to the long-range component :math:`C_{\mathrm{LR}}`. +CO\ :sub:`2` Concentration +===================================== + +The estimate of the concentration of CO\ :sub:`2` in a given room to indicate the air quality is given by the same approach as for the long-range virus concentration, +:math:`C_{\mathrm{LR}}(t, D)`, where :math:`C_0(D)` is considered to be the background (outdoor) CO\ :sub:`2` concentration (:meth:`caimira.models.CO2ConcentrationModel.CO2_atmosphere_concentration`). + +In order to compute the CO\ :sub:`2` concentration one should then simply use the :meth:`caimira.models.CO2ConcentrationModel.concentration` method. +A fraction of 4.2% of the exhalation rate of the defined activity was considered as supplied to the room (:meth:`caimira.models.CO2ConcentrationModel.CO2_fraction_exhaled`). + +Note still that nothing depends on the aerosol diameter :math:`D` in this case (no particles are involved) - hence in this class all parameters are constant w.r.t :math:`D`. + +Since the CO\ :sub:`2` concentration differs from the virus concentration, the specific removal rate, CO\ :sub:`2` atmospheric concentration and normalization factors are respectively defined in :meth:`caimira.models.CO2ConcentrationModel.removal_rate`, +:meth:`caimira.models.CO2ConcentrationModel.min_background_concentration` and :meth:`caimira.models.CO2ConcentrationModel.normalization_factor`. + .. _caimira-uml-diagram: CAiMIRA UML Diagram diff --git a/caimira/models.py b/caimira/models.py index 36da66f5..a9e5ae01 100644 --- a/caimira/models.py +++ b/caimira/models.py @@ -956,31 +956,43 @@ def probability_meet_infected_person(self, virus: Virus, n_infected: int, event_ @dataclass(frozen=True) -class ConcentrationModel: +class _ConcentrationModelBase: + """ + A generic superclass that contains the methods to calculate the + concentration (e.g. viral concentration or CO2 concentration). + """ room: Room ventilation: _VentilationBase - infected: InfectedPopulation - - #: evaporation factor: the particles' diameter is multiplied by this - # factor as soon as they are in the air (but AFTER going out of the, - # mask, if any). - evaporation_factor: float = 0.3 @property - def virus(self): - return self.infected.virus + def population(self) -> Population: + """ + Population in the room (the emitters of what we compute the + concentration of) + """ + raise NotImplementedError("Subclass must implement") - def infectious_virus_removal_rate(self, time: float) -> _VectorisedFloat: - # Equilibrium velocity of particle motion toward the floor - vg = self.infected.particle.settling_velocity(self.evaporation_factor) - # Height of the emission source to the floor - i.e. mouth/nose (m) - h = 1.5 - # Deposition rate (h^-1) - k = (vg * 3600) / h - return ( - k + self.virus.decay_constant(self.room.humidity, self.room.inside_temp.value(time)) - + self.ventilation.air_exchange(self.room, time) - ) + def removal_rate(self, time: float) -> _VectorisedFloat: + """ + Remove rate of the species considered, in h^-1 + """ + raise NotImplementedError("Subclass must implement") + + def min_background_concentration(self) -> _VectorisedFloat: + """ + Minimum background concentration in the room for a given scenario + (in the same unit as the concentration). Its the value towards which + the concentration will decay to. + """ + return 0. + + def normalization_factor(self) -> _VectorisedFloat: + """ + Normalization factor (in the same unit as the concentration). + This factor is applied to the normalized concentration only + at the very end. + """ + raise NotImplementedError("Subclass must implement") @method_cache def _normed_concentration_limit(self, time: float) -> _VectorisedFloat: @@ -988,27 +1000,27 @@ def _normed_concentration_limit(self, time: float) -> _VectorisedFloat: Provides a constant that represents the theoretical asymptotic value reached by the concentration when time goes to infinity, if all parameters were to stay time-independent. - This is normalized by the emission rate, the latter acting as a + This is normalized by the normalization factor, the latter acting as a multiplicative constant factor for the concentration model that can be put back in front of the concentration after the time dependence has been solved for. """ - if not self.infected.person_present(time): - return 0. + if not self.population.person_present(time): + return self.min_background_concentration()/self.normalization_factor() V = self.room.volume - IVRR = self.infectious_virus_removal_rate(time) + RR = self.removal_rate(time) - return 1. / (IVRR * V) + return (1. / (RR * V) + self.min_background_concentration()/ + self.normalization_factor()) @method_cache def state_change_times(self) -> typing.List[float]: """ All time dependent entities on this model must provide information about the times at which their state changes. - """ state_change_times = {0.} - state_change_times.update(self.infected.presence.transition_times()) + state_change_times.update(self.population.presence.transition_times()) state_change_times.update(self.ventilation.transition_times(self.room)) return sorted(state_change_times) @@ -1016,9 +1028,8 @@ def state_change_times(self) -> typing.List[float]: def _first_presence_time(self) -> float: """ First presence time. Before that, the concentration is zero. - """ - return self.infected.presence.boundaries()[0][0] + return self.population.presence.boundaries()[0][0] def last_state_change(self, time: float) -> float: """ @@ -1027,7 +1038,6 @@ def last_state_change(self, time: float) -> float: Find the nearest time less than the given one. If there is a state change exactly at ``time`` the previous state change is returned (except at ``time == 0``). - """ times = self.state_change_times() t_index: int = np.searchsorted(times, time) # type: ignore @@ -1040,7 +1050,6 @@ def last_state_change(self, time: float) -> float: def _next_state_change(self, time: float) -> float: """ Find the nearest future state change. - """ for change_time in self.state_change_times(): if change_time >= time: @@ -1052,15 +1061,17 @@ def _next_state_change(self, time: float) -> float: @method_cache def _normed_concentration_cached(self, time: float) -> _VectorisedFloat: - # A cached version of the _normed_concentration method. Use this - # method if you expect that there may be multiple concentration - # calculations for the same time (e.g. at state change times). + """ + A cached version of the _normed_concentration method. Use this + method if you expect that there may be multiple concentration + calculations for the same time (e.g. at state change times). + """ return self._normed_concentration(time) def _normed_concentration(self, time: float) -> _VectorisedFloat: """ - Virus long-range exposure concentration, as a function of time, and - normalized by the emission rate. + Concentration as a function of time, and normalized by + normalization_factor. The formulas used here assume that all parameters (ventilation, emission rate) are constant between two state changes - only the value of these parameters at the next state change, are used. @@ -1071,36 +1082,42 @@ def _normed_concentration(self, time: float) -> _VectorisedFloat: # 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 0.0 + return self.min_background_concentration()/self.normalization_factor() next_state_change_time = self._next_state_change(time) - IVRR = self.infectious_virus_removal_rate(next_state_change_time) - conc_limit = self._normed_concentration_limit(next_state_change_time) + RR = self.removal_rate(next_state_change_time) + # If RR is 0, conc_limit does not play a role but its computation + # would raise an error -> we set it to zero. + try: + conc_limit = self._normed_concentration_limit(next_state_change_time) + except ZeroDivisionError: + 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) delta_time = time - t_last_state_change - fac = np.exp(-IVRR * delta_time) + fac = np.exp(-RR * delta_time) return conc_limit * (1 - fac) + conc_at_last_state_change * fac def concentration(self, time: float) -> _VectorisedFloat: """ - Virus long-range exposure concentration, as a function of time. + Total concentration as a function of time. The normalization + factor has been put back. Note that time is not vectorised. You can only pass a single float to this method. """ return (self._normed_concentration_cached(time) * - self.infected.emission_rate_when_present()) + self.normalization_factor()) @method_cache def normed_integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat: """ - Get the integrated long-range concentration of viruses in the air between the times start and stop, - normalized by the emission rate. + Get the integrated concentration between the times start and stop, + normalized by normalization_factor. """ if stop <= self._first_presence_time(): - return 0.0 + return (stop - start)*self.min_background_concentration()/self.normalization_factor() state_change_times = self.state_change_times() req_start, req_stop = start, stop total_normed_concentration = 0. @@ -1115,11 +1132,11 @@ def normed_integrated_concentration(self, start: float, stop: float) -> _Vectori next_conc_state = self._next_state_change(stop) conc_limit = self._normed_concentration_limit(next_conc_state) - IVRR = self.infectious_virus_removal_rate(next_conc_state) + RR = self.removal_rate(next_conc_state) delta_time = stop - start total_normed_concentration += ( conc_limit * delta_time + - (conc_limit - conc_start) * (np.exp(-IVRR*delta_time)-1) / IVRR + (conc_limit - conc_start) * (np.exp(-RR*delta_time)-1) / RR ) return total_normed_concentration @@ -1128,7 +1145,85 @@ def integrated_concentration(self, start: float, stop: float) -> _VectorisedFloa Get the integrated concentration of viruses in the air between the times start and stop. """ return (self.normed_integrated_concentration(start, stop) * - self.infected.emission_rate_when_present()) + self.normalization_factor()) + + +@dataclass(frozen=True) +class ConcentrationModel(_ConcentrationModelBase): + """ + Class used for the computation of the long-range virus concentration. + """ + + #: Infected population in the room, emitting virions + infected: InfectedPopulation + + #: evaporation factor: the particles' diameter is multiplied by this + # factor as soon as they are in the air (but AFTER going out of the, + # mask, if any). + evaporation_factor: float = 0.3 + + @property + def population(self) -> InfectedPopulation: + return self.infected + + @property + def virus(self) -> Virus: + return self.infected.virus + + def normalization_factor(self) -> _VectorisedFloat: + # we normalize by the emission rate + return self.infected.emission_rate_when_present() + + def removal_rate(self, time: float) -> _VectorisedFloat: + # Equilibrium velocity of particle motion toward the floor + vg = self.infected.particle.settling_velocity(self.evaporation_factor) + # Height of the emission source to the floor - i.e. mouth/nose (m) + h = 1.5 + # Deposition rate (h^-1) + k = (vg * 3600) / h + return ( + k + self.virus.decay_constant(self.room.humidity, self.room.inside_temp.value(time)) + + self.ventilation.air_exchange(self.room, time) + ) + + def infectious_virus_removal_rate(self, time: float) -> _VectorisedFloat: + # defined for back-compatibility purposes + return self.removal_rate(time) + + +@dataclass(frozen=True) +class CO2ConcentrationModel(_ConcentrationModelBase): + """ + Class used for the computation of the CO2 concentration. + """ + #: Population in the room emitting CO2 + CO2_emitters: Population + + #: CO2 concentration in the atmosphere (in ppm) + CO2_atmosphere_concentration: float = 440.44 + + #: CO2 fraction in the exhaled air + CO2_fraction_exhaled: float = 0.042 + + @property + def population(self) -> Population: + return self.CO2_emitters + + def removal_rate(self, time: float) -> _VectorisedFloat: + return self.ventilation.air_exchange(self.room, time) + + def min_background_concentration(self) -> _VectorisedFloat: + """ + Background CO2 concentration in the atmosphere (in ppm) + """ + return self.CO2_atmosphere_concentration + + def normalization_factor(self) -> _VectorisedFloat: + # normalization by the CO2 exhaled. + # CO2 concentration given in ppm, hence the 1e6 factor. + return (1e6*self.population.number + *self.population.activity.exhalation_rate + *self.CO2_fraction_exhaled) @dataclass(frozen=True) @@ -1165,7 +1260,7 @@ def dilution_factor(self) -> _VectorisedFloat: φ = 2 # Exhalation airflow, as per Jia et al. (2022) - Q_exh = φ * BR + Q_exh: _VectorisedFloat = φ * BR # Area of the mouth assuming a perfect circle (m2) Am = np.pi*(mouth_diameter**2)/4 @@ -1349,14 +1444,13 @@ def __post_init__(self): """ c_model = self.concentration_model # Check if the diameter is vectorised. - if (isinstance(c_model.infected, InfectedPopulation) and not np.isscalar(c_model.infected.expiration.diameter) + if (isinstance(c_model.infected, InfectedPopulation) and not np.isscalar(c_model.infected.expiration.diameter) # Check if the diameter-independent elements of the infectious_virus_removal_rate method are vectorised. 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.") - def long_range_fraction_deposited(self) -> _VectorisedFloat: """ diff --git a/caimira/tests/models/test_co2_concentration_model.py b/caimira/tests/models/test_co2_concentration_model.py new file mode 100644 index 00000000..ad4a348c --- /dev/null +++ b/caimira/tests/models/test_co2_concentration_model.py @@ -0,0 +1,44 @@ +import numpy.testing as npt +import pytest + +from caimira import models + + +@pytest.fixture +def simple_co2_conc_model(): + return models.CO2ConcentrationModel( + room=models.Room(200, models.PiecewiseConstant((0., 24.), (293,))), + ventilation=models.AirChange(models.PeriodicInterval(period=120, duration=120), 0.25), + CO2_emitters=models.Population( + number=5, + presence=models.SpecificInterval((([0., 4.], ))), + mask=models.Mask.types['No mask'], + activity=models.Activity.types['Seated'], + host_immunity=0., + ), + ) + + +@pytest.mark.parametrize( + "time, expected_co2_concentration", [ + [0., 440.44], + [1., 914.2487227], + [2., 1283.251327], + [3., 1570.630844], + [4., 1794.442237], + ] +) +def test_co2_concentration( + simple_co2_conc_model: models.CO2ConcentrationModel, + time: float, + expected_co2_concentration: float, +): + npt.assert_almost_equal(simple_co2_conc_model.concentration(time), expected_co2_concentration) + + +def test_integrated_concentration(simple_co2_conc_model): + c1 = simple_co2_conc_model.integrated_concentration(0, 2) + c2 = simple_co2_conc_model.integrated_concentration(0, 1) + c3 = simple_co2_conc_model.integrated_concentration(1, 2) + assert c1 != 0 + npt.assert_almost_equal(c1, c2 + c3) diff --git a/caimira/tests/models/test_concentration_model.py b/caimira/tests/models/test_concentration_model.py index 0f6acb37..568e6757 100644 --- a/caimira/tests/models/test_concentration_model.py +++ b/caimira/tests/models/test_concentration_model.py @@ -3,9 +3,38 @@ import numpy as np import numpy.testing as npt import pytest +from dataclasses import dataclass from caimira import models +@dataclass(frozen=True) +class KnownConcentrationModelBase(models._ConcentrationModelBase): + """ + A _ConcentrationModelBase class where all the class methods are + redefined with a value taken from new parameters. Useful for testing. + + """ + known_population: models.Population + + known_removal_rate: float + + known_min_background_concentration: float + + known_normalization_factor: float + + @property + def population(self) -> models.Population: + return self.known_population + + def removal_rate(self, time: float) -> float: + return self.known_removal_rate + + def min_background_concentration(self) -> float: + return self.known_min_background_concentration + + def normalization_factor(self) -> float: + return self.known_normalization_factor + @pytest.mark.parametrize( "override_params", [ @@ -59,9 +88,9 @@ def test_concentration_model_vectorisation(override_params): def simple_conc_model(): interesting_times = models.SpecificInterval(([0.5, 1.], [1.1, 2], [2., 3.]), ) return models.ConcentrationModel( - models.Room(75, models.PiecewiseConstant((0., 24.), (293,))), - models.AirChange(interesting_times, 100), - models.InfectedPopulation( + room = models.Room(75, models.PiecewiseConstant((0., 24.), (293,))), + ventilation = models.AirChange(interesting_times, 100), + infected = models.InfectedPopulation( number=1, presence=interesting_times, mask=models.Mask.types['Type I'], @@ -74,6 +103,17 @@ def simple_conc_model(): ) +@pytest.fixture +def dummy_population(simple_conc_model) -> models.Population: + return models.Population( + number=10, + presence=simple_conc_model.infected.presence, + mask=models.Mask.types['Type I'], + activity=models.Activity.types['Seated'], + host_immunity=0., + ) + + @pytest.mark.parametrize( "time, expected_last_state_change", [ [-15., 0.], # Out of range goes to the first state. @@ -137,3 +177,98 @@ def test_integrated_concentration(simple_conc_model): c3 = simple_conc_model.integrated_concentration(1, 2) assert c1 != 0 npt.assert_almost_equal(c1, c2 + c3, decimal=15) + + +# The expected numbers were obtained via the quad integration of the +# normed_integrated_concentration method with 0 (start) and 2 (stop) as limits. +@pytest.mark.parametrize([ + "known_min_background_concentration", + "expected_normed_integrated_concentration"], + [ + [0.0, 0.00018533333708996207], + [240.0, 48.000185340695275], + [440.0, 88.00018534069527], + [600., 120.00018534069527], + [1000., 200.0001853407918], + ] +) +def test_normed_integrated_concentration_with_background_concentration( + simple_conc_model: models.ConcentrationModel, + dummy_population: models.Population, + known_min_background_concentration: float, + expected_normed_integrated_concentration: float): + + known_conc_model = KnownConcentrationModelBase( + room = simple_conc_model.room, + ventilation = simple_conc_model.ventilation, + known_population = dummy_population, + known_removal_rate = 100., + known_min_background_concentration = known_min_background_concentration, + known_normalization_factor = 10.) + npt.assert_almost_equal(known_conc_model.normed_integrated_concentration(0, 2), expected_normed_integrated_concentration) + + +# The expected numbers were obtained via the quad integration of the +# normed_integrated_concentration method with 0 (start) and 2 (stop) as limits. +@pytest.mark.parametrize([ + "known_removal_rate", + "known_min_background_concentration", + "known_normalization_factor", + "expected_normed_integrated_concentration"], + [ + [np.array([0.25, 10]), 0.0, 10., np.array([0.012161005755130391, 0.0017333437605308818])], + [100, np.array([0, 240.0]), 10., np.array([0.00018533333708996207, 48.000185340695275])], + [100, 440.0, np.array([10., 20.]), np.array([88.00018534069527, 44.000185340695275])], + [np.array([10, 100]), np.array([600., 800.]), 10., np.array([120.00173334473946, 160.0001853406953])], + [np.array([50, 100,]), np.array([1000.,1100.]), np.array([10., 20.]), np.array([200.00036800764332, 110.00018534069527])], + ] +) +def test_normed_integrated_concentration_vectorisation( + simple_conc_model: models.ConcentrationModel, + dummy_population: models.Population, + known_removal_rate: float, + known_min_background_concentration: float, + known_normalization_factor: float, + expected_normed_integrated_concentration: float): + + known_conc_model = KnownConcentrationModelBase( + room = simple_conc_model.room, + ventilation = simple_conc_model.ventilation, + known_population = dummy_population, + known_removal_rate = known_removal_rate, + known_min_background_concentration = known_min_background_concentration, + known_normalization_factor = known_normalization_factor) + + integrated_concentration = known_conc_model.normed_integrated_concentration(0, 2) + + assert isinstance(integrated_concentration, np.ndarray) + assert integrated_concentration.shape == (2, ) + npt.assert_almost_equal(integrated_concentration, expected_normed_integrated_concentration) + + +@pytest.mark.parametrize([ + "known_removal_rate", + "known_min_background_concentration", + "expected_concentration"], + [ + [0., 240., 240.], + [0., np.array([240., 240.]), np.array([240., 240.])] + ] +) +def test_zero_ventilation_rate( + simple_conc_model: models.ConcentrationModel, + dummy_population: models.Population, + known_removal_rate: float, + known_min_background_concentration: float, + expected_concentration: float): + + known_conc_model = KnownConcentrationModelBase( + room = simple_conc_model.room, + ventilation = simple_conc_model.ventilation, + known_population = dummy_population, + known_removal_rate = known_removal_rate, + known_normalization_factor=10., + known_min_background_concentration = known_min_background_concentration) + + normed_concentration = known_conc_model.concentration(1) + npt.assert_almost_equal(normed_concentration, expected_concentration) diff --git a/caimira/tests/models/test_exposure_model.py b/caimira/tests/models/test_exposure_model.py index 1a986815..ef9d0c4e 100644 --- a/caimira/tests/models/test_exposure_model.py +++ b/caimira/tests/models/test_exposure_model.py @@ -19,7 +19,7 @@ class KnownNormedconcentration(models.ConcentrationModel): """ normed_concentration_function: typing.Callable = lambda x: 0 - def infectious_virus_removal_rate(self, time: float) -> models._VectorisedFloat: + def removal_rate(self, time: float) -> models._VectorisedFloat: # Very large decay constant -> same as constant concentration return 1.e50 diff --git a/caimira/tests/test_full_algorithm.py b/caimira/tests/test_full_algorithm.py index 51a3da72..9f2ee70d 100644 --- a/caimira/tests/test_full_algorithm.py +++ b/caimira/tests/test_full_algorithm.py @@ -229,7 +229,7 @@ def dilution_factor(self) -> _VectorisedFloat: x = np.array(self.distance) dilution = np.empty(x.shape, dtype=np.float64) # Exhalation airflow, as per Jia et al. (2022), m^3/s - Q_exh = self.φ * np.array(self.breathing_rate/3600) + Q_exh: _VectorisedFloat = self.φ * np.array(self.breathing_rate/3600) # The expired flow velocity at the noozle (mouth opening), m/s u0 = np.array(Q_exh/(np.pi/4. * self.mouth_diameter**2)) # Parameters in the jet-like stage