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

Merge branch 'feature/concentration_normalization' into 'master'

Browse source 
Normalization of the computation of the concentration

See merge request cara/cara!255
This commit is contained in:
Andre Henriques 2021-09-14 16:11:45 +02:00
commit 705485b1d5
3 changed files with 70 additions and 42 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

94
cara/models.py
View file

@ -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()

16
cara/tests/models/test_exposure_model.py
View file

@ -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(

2
cara/tests/test_monte_carlo_full_models.py
View file

@ -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],
]