From 43da2d7521479ce0609c3bf60f3cebec9179cf89 Mon Sep 17 00:00:00 2001 From: Phil Elson Date: Tue, 10 Nov 2020 15:45:39 +0100 Subject: [PATCH] Implement exposed activity/mask wearing independently of the infected group. --- cara/models.py | 84 +++++++++++++++++++++-------- cara/tests/test_known_quantities.py | 82 ++++++++++++++++++---------- 2 files changed, 116 insertions(+), 50 deletions(-) diff --git a/cara/models.py b/cara/models.py index 61a2182d..d765b548 100644 --- a/cara/models.py +++ b/cara/models.py @@ -389,19 +389,41 @@ Activity.types = { @dataclass(frozen=True) -class InfectedPerson: - virus: Virus - #: The times in which the person is in the room. +class Population: + """ + Represents a group of people all with exactly the same behaviour and + situation. + + """ + #: How many in the population. + number: int + + #: The times in which the people are in the room. presence: Interval + + #: The kind of mask being worn by the people. mask: Mask + + #: The physical activity being carried out by the people. activity: Activity - expiration: Expiration def person_present(self, time): return self.presence.triggered(time) - @functools.lru_cache() - def emission_rate(self, time) -> float: + +@dataclass(frozen=True) +class InfectedPopulation(Population): + #: The virus with which the population is infected. + virus: Virus + + #: The type of expiration that is being emitted whilst doing the activity. + expiration: Expiration + + def individual_emission_rate(self, time) -> float: + """ + The emission rate of a single individual in the population. + + """ # Note: The original model avoids time dependence on the emission rate # at the cost of implementing a piecewise (on time) concentration function. if not self.person_present(time): @@ -420,15 +442,20 @@ class InfectedPerson: aerosols) return ER + @functools.lru_cache() + def emission_rate(self, time) -> float: + """ + The emission rate of the entire population. + + """ + return self.individual_emission_rate(time) * self.number + @dataclass(frozen=True) class Model: room: Room ventilation: Ventilation - infected: InfectedPerson - infected_occupants: int - exposed_occupants: int - exposed_activity: Activity + infected: InfectedPopulation @property def virus(self): @@ -473,36 +500,49 @@ class Model: return 0.0 IVRR = self.infectious_virus_removal_rate(time) V = self.room.volume - Ni = self.infected_occupants - ER = self.infected.emission_rate(time) t_last_state_change = self.last_state_change(time) concentration_at_last_state_change = self.concentration(t_last_state_change) delta_time = time - t_last_state_change fac = np.exp(-IVRR * delta_time) - concentration_limit = (ER * Ni) / (IVRR * V) + concentration_limit = (self.infected.emission_rate(time)) / (IVRR * V) return concentration_limit * (1 - fac) + concentration_at_last_state_change * fac - def infection_probability(self): - # Infection probability - # Probability of COVID-19 Infection - exposure = 0.0 # q/m3*h +@dataclass(frozen=True) +class ExposureModel: + #: The virus concentration model which this exposure model should consider. + concentration_model: Model + + #: The population of non-infected people to be used in the model. + exposed: Population + + def quanta_exposure(self) -> float: + """The number of virus quanta per meter^3.""" + exposure = 0.0 def integrate(fn, start, stop): values = np.linspace(start, stop) return np.trapz([fn(v) for v in values], values) - # TODO: Have this for exposed not infected. - for start, stop in self.infected.presence.boundaries(): - exposure += (integrate(self.concentration, start, stop)) + for start, stop in self.exposed.presence.boundaries(): + exposure += integrate(self.concentration_model.concentration, start, stop) + return exposure + + def infection_probability(self): + exposure = self.quanta_exposure() inf_aero = ( - self.exposed_activity.inhalation_rate * - (1 - self.infected.mask.η_inhale) * + self.exposed.activity.inhalation_rate * + (1 - self.exposed.mask.η_inhale) * exposure ) # Probability of infection. return (1 - np.exp(-inf_aero)) * 100 + + def reproduction_rate(self): + prob = self.infection_probability() + exposed_occupants = self.exposed.number + return prob * exposed_occupants / 100 diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 29d16d72..30e64447 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -34,20 +34,31 @@ def baseline_model(): outside_temp=models.PiecewiseConstant((0,24),(283,)), cd_b=0.6, window_height=1.6, opening_length=0.6, ), - infected=models.InfectedPerson( + infected=models.InfectedPopulation( + number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(((0, 4), (5, 8))), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light exercise'], expiration=models.Expiration.types['Unmodulated Vocalization'], ), - infected_occupants=1, - exposed_occupants=10, - exposed_activity=models.Activity.types['Light exercise'], ) return model +@pytest.fixture +def baseline_exposure_model(baseline_model): + return models.ExposureModel( + baseline_model, + exposed=models.Population( + number=10, + presence=baseline_model.infected.presence, + activity=baseline_model.infected.activity, + mask=baseline_model.infected.mask, + ) + ) + + @pytest.fixture def baseline_periodic_window(): return models.WindowOpening( @@ -98,16 +109,14 @@ def build_model(interval_duration): active=models.PeriodicInterval(period=120, duration=interval_duration), q_air_mech=500., ), - infected=models.InfectedPerson( + infected=models.InfectedPopulation( + number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(((0, 4), (5, 8))), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light exercise'], expiration=models.Expiration.types['Unmodulated Vocalization'], ), - infected_occupants=1, - exposed_occupants=10, - exposed_activity=models.Activity.types['Light exercise'], ) return model @@ -121,8 +130,8 @@ def test_concentrations_startup(baseline_model): assert m1.concentration(1.) == m2.concentration(1.) -def test_r0(baseline_model): - p = baseline_model.infection_probability() +def test_r0(baseline_exposure_model): + p = baseline_exposure_model.infection_probability() npt.assert_allclose(p, 93.196908) @@ -170,7 +179,6 @@ def test_multiple_ventilation_HEPA_window(baseline_periodic_hepa, time, expected def test_multiple_ventilation_HEPA_window_transitions(baseline_periodic_hepa): - room = models.Room(volume=68.) tempOutside = models.PiecewiseConstant((0., 1., 2.5),(273.15, 283.15)) tempInside = models.PiecewiseConstant((0., 24.),(293.15,)) window = models.WindowOpening(active=models.SpecificInterval([(1 / 60, 24.)]), @@ -322,16 +330,14 @@ def build_hourly_dependent_model(month, intervals_open=((7.5, 8.5),), outside_temp=outside_temp, cd_b=0.6, window_height=1.6, opening_length=0.6, ), - infected=models.InfectedPerson( + infected=models.InfectedPopulation( + number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(intervals_presence_infected), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light exercise'], expiration=models.Expiration.types['Unmodulated Vocalization'], ), - infected_occupants=1, - exposed_occupants=10, - exposed_activity=models.Activity.types['Light exercise'], ) return model @@ -345,16 +351,14 @@ def build_constant_temp_model(outside_temp, intervals_open=((7.5, 8.5),)): outside_temp=models.PiecewiseConstant((0,24),(outside_temp,)), cd_b=0.6, window_height=1.6, opening_length=0.6, ), - infected=models.InfectedPerson( + infected=models.InfectedPopulation( + number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(((0, 4), (5, 7.5))), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light exercise'], expiration=models.Expiration.types['Unmodulated Vocalization'], ), - infected_occupants=1, - exposed_occupants=10, - exposed_activity=models.Activity.types['Light exercise'], ) return model @@ -374,16 +378,14 @@ def build_hourly_dependent_model_multipleventilation(month, intervals_open=((7.5 model = models.Model( room=models.Room(volume=75), ventilation=vent, - infected=models.InfectedPerson( + infected=models.InfectedPopulation( + number=1, virus=models.Virus.types['SARS_CoV_2'], presence=models.SpecificInterval(((0, 4), (5, 7.5))), mask=models.Mask.types['No mask'], activity=models.Activity.types['Light exercise'], expiration=models.Expiration.types['Unmodulated Vocalization'], ), - infected_occupants=1, - exposed_occupants=10, - exposed_activity=models.Activity.types['Light exercise'], ) return model @@ -451,6 +453,20 @@ def test_concentrations_refine_times(time): artificial_refinement=True) npt.assert_allclose(m1.concentration(time), m2.concentration(time), rtol=1e-8) + +def build_exposure_model(concentration_model): + infected = concentration_model.infected + return models.ExposureModel( + concentration_model=concentration_model, + exposed=models.Population( + number=10, + presence=infected.presence, + activity=infected.activity, + mask=infected.mask, + ) + ) + + @pytest.mark.parametrize( "month, expected_r0", [ @@ -459,8 +475,13 @@ def test_concentrations_refine_times(time): ], ) def test_r0_hourly_dep(month,expected_r0): - m = build_hourly_dependent_model(month,intervals_open=((0,24),), - intervals_presence_infected=((8,12),(13,17))) + m = build_exposure_model( + build_hourly_dependent_model( + month, + intervals_open=((0,24),), + intervals_presence_infected=((8, 12), (13, 17)) + ) + ) p = m.infection_probability() npt.assert_allclose(p, expected_r0) @@ -472,8 +493,13 @@ def test_r0_hourly_dep(month,expected_r0): ], ) def test_r0_hourly_dep_refined(month,expected_r0): - m = build_hourly_dependent_model(month,intervals_open=((0,24),), - intervals_presence_infected=((8,12),(13,17)), - temperatures=data.GenevaTemperatures) + m = build_exposure_model( + build_hourly_dependent_model( + month, + intervals_open=((0, 24),), + intervals_presence_infected=((8, 12), (13, 17)), + temperatures=data.GenevaTemperatures, + ) + ) p = m.infection_probability() npt.assert_allclose(p, expected_r0)