diff --git a/cara/apps/calculator/model_generator.py b/cara/apps/calculator/model_generator.py index 7e4a95c3..27f10ad5 100644 --- a/cara/apps/calculator/model_generator.py +++ b/cara/apps/calculator/model_generator.py @@ -120,7 +120,7 @@ class FormData: infected_finish=time_string_to_minutes(form_data['infected_finish']), ) - def build_model(self) -> models.Model: + def build_model(self) -> models.ExposureModel: return model_from_form(self) def ventilation(self) -> models.Ventilation: @@ -224,7 +224,7 @@ class FormData: return models.SpecificInterval(tuple(present_intervals)) -def model_from_form(form: FormData) -> models.Model: +def model_from_form(form: FormData) -> models.ExposureModel: # Initializes room with volume either given directly or as product of area and height if form.volume_type == 'room_volume': volume = form.room_volume @@ -257,19 +257,25 @@ def model_from_form(form: FormData) -> models.Model: exposed_occupants = form.total_people - infected_occupants # Initializes and returns a model with the attributes defined above - return models.Model( - room=room, - ventilation=form.ventilation(), - infected=models.InfectedPerson( - virus=virus, - presence=form.present_interval(), - mask=mask, - activity=infected_activity, - expiration=infected_expiration + return models.ExposureModel( + concentration_model=models.ConcentrationModel( + room=room, + ventilation=form.ventilation(), + infected=models.InfectedPopulation( + number=infected_occupants, + virus=virus, + presence=form.present_interval(), + mask=mask, + activity=infected_activity, + expiration=infected_expiration + ), ), - infected_occupants=infected_occupants, - exposed_occupants=exposed_occupants, - exposed_activity=exposed_activity + exposed=models.Population( + number=exposed_occupants, + presence=form.present_interval(), + activity=exposed_activity, + mask=mask, + ) ) diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index ee9aa22e..2a8fa1c5 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -13,20 +13,19 @@ from cara import models from .model_generator import FormData -def calculate_report_data(model: models.Model): +def calculate_report_data(model: models.ExposureModel): resolution = 600 - # TODO: Have this for exposed not infected. - t_start = model.infected.presence.boundaries()[0][0] - t_end = model.infected.presence.boundaries()[-1][1] + t_start = model.exposed.presence.boundaries()[0][0] + t_end = model.exposed.presence.boundaries()[-1][1] times = list(np.linspace(t_start, t_end, resolution)) - concentrations = [model.concentration(time) for time in times] + concentrations = [model.concentration_model.concentration(time) for time in times] highest_const = max(concentrations) prob = model.infection_probability() - er = model.infected.emission_rate(0.1) - exposed_occupants = model.exposed_occupants - r0 = prob * exposed_occupants / 100 + er = model.concentration_model.infected.emission_rate_when_present() + exposed_occupants = model.exposed.number + r0 = model.reproduction_rate() return { "times": times, @@ -78,7 +77,7 @@ def minutes_to_time(minutes: int) -> str: return f"{hour_string}:{minute_string}" -def build_report(model: models.Model, form: FormData): +def build_report(model: models.ExposureModel, form: FormData): now = datetime.now() time = now.strftime("%d/%m/%Y %H:%M:%S") request = {"the": "form", "request": "data"} diff --git a/cara/apps/calculator/templates/report.html.j2 b/cara/apps/calculator/templates/report.html.j2 index 4a4ccfda..0be69661 100644 --- a/cara/apps/calculator/templates/report.html.j2 +++ b/cara/apps/calculator/templates/report.html.j2 @@ -20,7 +20,7 @@

Input data:

Ventilation data:

diff --git a/cara/apps/expert.py b/cara/apps/expert.py index c02c56df..0d8009ba 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -33,7 +33,7 @@ class ConcentrationFigure: self.ax = self.figure.add_subplot(1, 1, 1) self.line = None - def update(self, model: models.Model): + def update(self, model: models.ConcentrationModel): resolution = 600 ts = np.linspace(0, 10, resolution) concentration = [model.concentration(t) for t in ts] @@ -105,30 +105,29 @@ class WidgetView: pass def update(self): - model = self.model_state.dcs_instance() + model: models.ExposureModel = self.model_state.dcs_instance() for plot in self.plots: - plot.update(model) + plot.update(model.concentration_model) self.out.clear_output() with self.out: P = model.infection_probability() - print(f'Emission rate (quanta/hr): {model.infected.emission_rate(0)}') + print(f'Emission rate (quanta/hr): {model.concentration_model.infected.emission_rate_when_present()}') print(f'Probability of infection: {np.round(P, 0)}%') - print(f'Number of exposed: {model.exposed_occupants}') - R0 = np.round(P / 100 * model.exposed_occupants, 1) + print(f'Number of exposed: {model.exposed.number}') + R0 = np.round(model.reproduction_rate(), 1) print(f'Number of expected new cases (R0): {R0}') - def _build_widget(self, node): - self.widget.children += (self._build_room(node.room),) - self.widget.children += (self._build_ventilation(node.ventilation),) - self.widget.children += (self._build_infected(node.infected),) + self.widget.children += (self._build_room(node.concentration_model.room),) + self.widget.children += (self._build_ventilation(node.concentration_model.ventilation),) + self.widget.children += (self._build_infected(node.concentration_model.infected),) self.widget.children += (self._build_exposed(node),) def _build_exposed(self, node): return collapsible( - [self._build_activity(node.exposed_activity)], + [self._build_activity(node.exposed.activity)], title="Exposed" ) @@ -284,24 +283,30 @@ class WidgetView: return self.widget -baseline_model = models.Model( - room=models.Room(volume=75), - ventilation=models.WindowOpening( - active=models.PeriodicInterval(period=120, duration=120), - inside_temp=models.PiecewiseConstant((0,24),(293,)), - outside_temp=models.PiecewiseConstant((0,24),(283,)), - cd_b=0.6, window_height=1.6, opening_length=0.6, +baseline_model = models.ExposureModel( + concentration_model=models.ConcentrationModel( + room=models.Room(volume=75), + ventilation=models.WindowOpening( + active=models.PeriodicInterval(period=120, duration=120), + inside_temp=models.PiecewiseConstant((0,24),(293,)), + outside_temp=models.PiecewiseConstant((0,24),(283,)), + cd_b=0.6, window_height=1.6, opening_length=0.6, + ), + 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=models.InfectedPerson( - virus=models.Virus.types['SARS_CoV_2'], + exposed=models.Population( + number=10, 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'], + mask=models.Mask.types['No mask'], ), - infected_occupants=1, - exposed_occupants=10, - exposed_activity=models.Activity.types['Light exercise'], ) @@ -330,13 +335,13 @@ class CARAStateBuilder(state.StateBuilder): class ExpertApplication: def __init__(self): self.model_state = state.DataclassInstanceState( - models.Model, + models.ExposureModel, state_builder=CARAStateBuilder(), ) self.model_state.dcs_update_from(baseline_model) # For the time-being, we have to initialise the select states. Careful # as values might not correspond to what the baseline model says. - self.model_state.infected.mask.dcs_select('No mask') + self.model_state.concentration_model.infected.mask.dcs_select('No mask') self.view = WidgetView(self.model_state) diff --git a/cara/models.py b/cara/models.py index 10ff25c8..2fcd576d 100644 --- a/cara/models.py +++ b/cara/models.py @@ -394,24 +394,43 @@ 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: - # 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): - return 0 +@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 emission_rate_when_present(self) -> float: + """ + The emission rate if the infected population is present. + + Note that the rate is not currently time-dependent. + + """ # Emission Rate (infectious quantum / h) aerosols = self.expiration.aerosols(self.mask) if np.isinf(aerosols): @@ -425,15 +444,38 @@ class InfectedPerson: aerosols) return ER + 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): + return 0. + + # Note: It is essential that the value of the emission rate is not + # 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() + + @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: +class ConcentrationModel: room: Room ventilation: Ventilation - infected: InfectedPerson - infected_occupants: int - exposed_occupants: int - exposed_activity: Activity + infected: InfectedPopulation @property def virus(self): @@ -478,36 +520,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: ConcentrationModel + + #: 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/apps/test_expert_app.py b/cara/tests/apps/test_expert_app.py index 26851489..af441e99 100644 --- a/cara/tests/apps/test_expert_app.py +++ b/cara/tests/apps/test_expert_app.py @@ -6,4 +6,4 @@ def test_app(): # do anything fancy to verify how it looks etc., we leave that for manual # testing. expert_app = cara.apps.ExpertApplication() - assert expert_app.model_state.room.volume == 75 + assert expert_app.model_state.concentration_model.room.volume == 75 diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 29d16d72..336751f9 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -26,7 +26,7 @@ def test_no_mask_emission_rate(baseline_model): @pytest.fixture def baseline_model(): - model = models.Model( + model = models.ConcentrationModel( room=models.Room(volume=75), ventilation=models.WindowOpening( active=models.PeriodicInterval(period=120, duration=120), @@ -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( @@ -92,22 +103,20 @@ def test_smooth_concentrations(baseline_model): def build_model(interval_duration): - model = models.Model( + model = models.ConcentrationModel( room=models.Room(volume=75), ventilation=models.HEPAFilter( 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.)]), @@ -314,7 +322,7 @@ def build_hourly_dependent_model(month, intervals_open=((7.5, 8.5),), else: outside_temp = temperatures[month] - model = models.Model( + model = models.ConcentrationModel( room=models.Room(volume=75), ventilation=models.WindowOpening( active=models.SpecificInterval(intervals_open), @@ -322,22 +330,20 @@ 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 def build_constant_temp_model(outside_temp, intervals_open=((7.5, 8.5),)): - model = models.Model( + model = models.ConcentrationModel( room=models.Room(volume=75), ventilation=models.WindowOpening( active=models.SpecificInterval(intervals_open), @@ -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 @@ -371,19 +375,17 @@ def build_hourly_dependent_model_multipleventilation(month, intervals_open=((7.5 active=models.SpecificInterval(((0,24),)), q_air_mech=500., ))) - model = models.Model( + model = models.ConcentrationModel( 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)