diff --git a/cara/models.py b/cara/models.py index 86d2acb9..98245997 100644 --- a/cara/models.py +++ b/cara/models.py @@ -701,11 +701,11 @@ class InfectedPopulation(Population): elif np.isinf(aerosols): ER = 970 * self.virus.infectious_dose - return ER + return ER * self.number - def individual_emission_rate(self, time) -> _VectorisedFloat: + def emission_rate(self, time) -> _VectorisedFloat: """ - The emission rate of a single individual in the population. + The emission rate of the population. """ # Note: The original model avoids time dependence on the emission rate @@ -721,13 +721,6 @@ class InfectedPopulation(Population): return self.emission_rate_when_present() - def emission_rate(self, time) -> _VectorisedFloat: - """ - The emission rate of the entire population. - - """ - return self.individual_emission_rate(time) * self.number - @dataclass(frozen=True) class ConcentrationModel: @@ -753,16 +746,22 @@ class ConcentrationModel: ) @method_cache - def _concentration_limit(self, time: float) -> _VectorisedFloat: + 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 + 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. V = self.room.volume IVRR = self.infectious_virus_removal_rate(time) - return (self.infected.emission_rate(time)) / (IVRR * V) + return 1. / (IVRR * V) @method_cache def state_change_times(self) -> typing.List[float]: @@ -815,15 +814,16 @@ class ConcentrationModel: ) @method_cache - def _concentration_cached(self, time: float) -> _VectorisedFloat: - # A cached version of the 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.concentration(time) + 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). + return self._normed_concentration(time) - def concentration(self, time: float) -> _VectorisedFloat: + def _normed_concentration(self, time: float) -> _VectorisedFloat: """ - Virus exposure concentration, as a function of time. + Virus exposure concentration, as a function of time, and + normalized by the emission rate. 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. @@ -837,25 +837,36 @@ class ConcentrationModel: return 0.0 next_state_change_time = self._next_state_change(time) IVRR = self.infectious_virus_removal_rate(next_state_change_time) - concentration_limit = self._concentration_limit(next_state_change_time) + conc_limit = self._normed_concentration_limit(next_state_change_time) t_last_state_change = self.last_state_change(time) - concentration_at_last_state_change = self._concentration_cached(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(-IVRR * delta_time) - return concentration_limit * (1 - fac) + concentration_at_last_state_change * fac + return conc_limit * (1 - fac) + conc_at_last_state_change * fac + + def concentration(self, time: float) -> _VectorisedFloat: + """ + Virus exposure concentration, as a function of time. + + Note that time is not vectorised. You can only pass a single float + to this method. + """ + return (self._normed_concentration(time) * + self.infected.emission_rate_when_present()) @method_cache - def integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat: + def normed_integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat: """ - Get the integrated concentration dose between the times start and stop. + Get the integrated concentration dose between the times start and stop, + normalized by the emission rate. """ if stop <= self._first_presence_time(): return 0.0 state_change_times = self.state_change_times() req_start, req_stop = start, stop - total_concentration = 0. + total_normed_concentration = 0. for interval_start, interval_stop in zip(state_change_times[:-1], state_change_times[1:]): if req_start > interval_stop or req_stop < interval_start: continue @@ -863,17 +874,24 @@ class ConcentrationModel: start = max([interval_start, req_start]) stop = min([interval_stop, req_stop]) - conc_start = self._concentration_cached(start) + conc_start = self._normed_concentration_cached(start) next_conc_state = self._next_state_change(stop) - conc_limit = self._concentration_limit(next_conc_state) + conc_limit = self._normed_concentration_limit(next_conc_state) IVRR = self.infectious_virus_removal_rate(next_conc_state) delta_time = stop - start - total_concentration += ( + total_normed_concentration += ( conc_limit * delta_time + (conc_limit - conc_start) * (np.exp(-IVRR*delta_time)-1) / IVRR ) - return total_concentration + return total_normed_concentration + + def integrated_concentration(self, start: float, stop: float) -> _VectorisedFloat: + """ + Get the integrated concentration dose between the times start and stop. + """ + return (self.normed_integrated_concentration(start, stop) * + self.infected.emission_rate_when_present()) @dataclass(frozen=True) @@ -890,14 +908,22 @@ class ExposureModel: #: The fraction of viruses actually deposited in the respiratory tract fraction_deposited: _VectorisedFloat = 0.6 - def exposure(self) -> _VectorisedFloat: - """The number of virus per meter^3.""" - exposure = 0.0 + def _normed_exposure(self) -> _VectorisedFloat: + """ + The number of virus per meter^3, normalized by the emission rate + of the infected population. + """ + normed_exposure = 0.0 for start, stop in self.exposed.presence.boundaries(): - exposure += self.concentration_model.integrated_concentration(start, stop) + normed_exposure += self.concentration_model.normed_integrated_concentration(start, stop) - return exposure * self.repeats + return normed_exposure * self.repeats + + def exposure(self) -> _VectorisedFloat: + """The number of virus per meter^3.""" + return (self._normed_exposure() * + self.concentration_model.infected.emission_rate_when_present()) def infection_probability(self) -> _VectorisedFloat: exposure = self.exposure() diff --git a/cara/tests/models/test_exposure_model.py b/cara/tests/models/test_exposure_model.py index cf8cfb6e..11a135a1 100644 --- a/cara/tests/models/test_exposure_model.py +++ b/cara/tests/models/test_exposure_model.py @@ -11,20 +11,20 @@ from cara.dataclass_utils import replace @dataclass(frozen=True) -class KnownConcentrations(models.ConcentrationModel): +class KnownNormedconcentration(models.ConcentrationModel): """ A ConcentrationModel which is based on pre-known exposure concentrations and which therefore doesn't need other components. Useful for testing. """ - concentration_function: typing.Callable + normed_concentration_function: typing.Callable def infectious_virus_removal_rate(self, time: float) -> models._VectorisedFloat: # very large decay constant -> same as constant concentration return 1.e50 - def _concentration_limit(self, time: float) -> models._VectorisedFloat: - return self.concentration_function(time) + def _normed_concentration_limit(self, time: float) -> models._VectorisedFloat: + return self.normed_concentration_function(time) def state_change_times(self): return [0., 24.] @@ -32,8 +32,8 @@ class KnownConcentrations(models.ConcentrationModel): def _next_state_change(self, time: float): return 24. - def concentration(self, time: float) -> models._VectorisedFloat: # noqa - return self.concentration_function(time) + def _normed_concentration(self, time: float) -> models._VectorisedFloat: # noqa + return self.normed_concentration_function(time) halftime = models.PeriodicInterval(120, 60) @@ -67,7 +67,9 @@ def known_concentrations(func): virus=models.Virus.types['SARS_CoV_2_B117'], expiration=models.Expiration.types['Talking'] ) - return KnownConcentrations(dummy_room, dummy_ventilation, dummy_infected_population, func) + normed_func = lambda x: func(x) / dummy_infected_population.emission_rate_when_present() + return KnownNormedconcentration(dummy_room, dummy_ventilation, + dummy_infected_population, normed_func) @pytest.mark.parametrize( diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index 5cf88198..08bd2dad 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -238,7 +238,7 @@ def skagit_chorale_mc(): ["shared_office_mc", 10.7, 0.32, 57.24, 654], ["classroom_mc", 36.1, 6.85, 780.0, 28464], ["ski_cabin_mc", 16.3, 0.49, 35.94, 7404], - ["gym_mc", 2.25, 0.63, 0.7842, 984], + ["gym_mc", 2.25, 0.63, 0.7842, 1968], ["waiting_room_mc", 9.72, 1.36, 34.26, 3534], ["skagit_chorale_mc",29.9, 17.9, 190.0, 141400], ]