From 31f628f9e5bce2c51382895e9ccc5b5e6e10c480 Mon Sep 17 00:00:00 2001 From: Phil Elson Date: Thu, 12 Nov 2020 21:21:02 +0100 Subject: [PATCH] Implement an ExposureModel.reproduction_number, which is modelled by having precisely one infected person in the concentration model. --- cara/apps/calculator/report_generator.py | 8 ++-- cara/apps/calculator/templates/report.html.j2 | 5 ++- cara/apps/expert.py | 8 +++- cara/dataclass_utils.py | 30 +++++++++++++++ cara/models.py | 26 ++++++++++--- cara/tests/conftest.py | 38 +++++++++++++++++++ cara/tests/test_dataclass_utils.py | 27 +++++++++++++ cara/tests/test_known_quantities.py | 35 ----------------- cara/tests/test_model.py | 12 ++++++ 9 files changed, 140 insertions(+), 49 deletions(-) create mode 100644 cara/dataclass_utils.py create mode 100644 cara/tests/conftest.py create mode 100644 cara/tests/test_dataclass_utils.py create mode 100644 cara/tests/test_model.py diff --git a/cara/apps/calculator/report_generator.py b/cara/apps/calculator/report_generator.py index b78d5069..0f05e1e9 100644 --- a/cara/apps/calculator/report_generator.py +++ b/cara/apps/calculator/report_generator.py @@ -18,7 +18,7 @@ from .model_generator import FormData class RepeatEvents: repeats: int probability_of_infection: float - R0: float + expected_new_cases: float def calculate_report_data(model: models.ExposureModel): @@ -33,7 +33,7 @@ def calculate_report_data(model: models.ExposureModel): prob = model.infection_probability() er = model.concentration_model.infected.emission_rate_when_present() exposed_occupants = model.exposed.number - r0 = model.reproduction_rate() + expected_new_cases = model.expected_new_cases() repeated_events = [] for n in [1, 2, 3, 4, 5, 10, 15, 20]: @@ -42,7 +42,7 @@ def calculate_report_data(model: models.ExposureModel): RepeatEvents( repeats=n, probability_of_infection=repeat_model.infection_probability(), - R0=repeat_model.reproduction_rate(), + expected_new_cases=repeat_model.expected_new_cases(), ) ) @@ -53,7 +53,7 @@ def calculate_report_data(model: models.ExposureModel): "prob_inf": prob, "emission_rate": er, "exposed_occupants": exposed_occupants, - "R0": r0, + "expected_new_cases": expected_new_cases, "scenario_plot_src": embed_figure(plot(times, concentrations)), "repeated_events": repeated_events, } diff --git a/cara/apps/calculator/templates/report.html.j2 b/cara/apps/calculator/templates/report.html.j2 index be04900f..c24b7a95 100644 --- a/cara/apps/calculator/templates/report.html.j2 +++ b/cara/apps/calculator/templates/report.html.j2 @@ -124,7 +124,8 @@

Results:

- In this scenario, the estimated probability of one exposed occupant getting infected P(i) is {{ prob_inf | int_format }}% and the expected number of new cases is {{ R0 | float_format }}. + In this scenario, the estimated probability of one exposed occupant getting infected P(i) is {{ prob_inf | int_format }}% + and the expected number of new cases is {{ expected_new_cases | float_format }}.

Exposure graph:

@@ -146,7 +147,7 @@ {{ repeat_event.repeats }} {{ repeat_event.probability_of_infection | int_format }}% - {{ repeat_event.R0 | float_format }} + {{ repeat_event.expected_new_cases | float_format }} {% endfor %} diff --git a/cara/apps/expert.py b/cara/apps/expert.py index 3a5f0326..60fcfee9 100644 --- a/cara/apps/expert.py +++ b/cara/apps/expert.py @@ -118,8 +118,12 @@ class WidgetView: print(f'Probability of infection: {np.round(P, 0)}%') print(f'Number of exposed: {model.exposed.number}') - R0 = np.round(model.reproduction_rate(), 1) - print(f'Number of expected new cases (R0): {R0}') + + new_cases = np.round(model.expected_new_cases(), 1) + print(f'Number of expected new cases: {new_cases}') + + R0 = np.round(model.reproduction_number(), 1) + print(f'Reproduction number (R0): {R0}') def _build_widget(self, node): self.widget.children += (self._build_room(node.concentration_model.room),) diff --git a/cara/dataclass_utils.py b/cara/dataclass_utils.py new file mode 100644 index 00000000..50654d43 --- /dev/null +++ b/cara/dataclass_utils.py @@ -0,0 +1,30 @@ +import dataclasses +import typing + + +DCInst = typing.TypeVar('T') + + +def nested_replace(obj: DCInst, new_values: typing.Dict[str, typing.Any]) -> DCInst: + """ + Replace an attribute on a dataclass, much like dataclasses.replace, except it + supports nested replacement definitions. For example: + + >>> new_obj = nested_replace(obj, {'attr1.sub_attr2.sub_sub_attr3': 4}) + >>> new_obj.attr1.sub_attr2.sub_sub_attr3 + 4 + + """ + new_inst = obj + for name, value in new_values.items(): + if '.' in name: + # Recurse into the desired name and come out with a top-level + # dataclass which has been updated appropriately. + name, remainder = name.split('.', 1) + value = nested_replace( + getattr(new_inst, name), + {remainder: value} + ) + # We have a plain old name. So set it. + new_inst = dataclasses.replace(new_inst, **{name: value}) + return new_inst diff --git a/cara/models.py b/cara/models.py index 12f85c7c..3825ef4b 100644 --- a/cara/models.py +++ b/cara/models.py @@ -1,10 +1,9 @@ +from dataclasses import dataclass import functools import numpy as np import typing -from abc import abstractmethod - -from dataclasses import dataclass +from .dataclass_utils import nested_replace @dataclass(frozen=True) @@ -138,7 +137,6 @@ class Ventilation: def transition_times(self) -> typing.Set[float]: return self.active.transition_times() - @abstractmethod def air_exchange(self, room: Room, time: float) -> float: """ Returns the rate at which air is being exchanged in the given room @@ -169,7 +167,6 @@ class MultipleVentilation: transitions.update(ventilation.transition_times()) return transitions - @abstractmethod def air_exchange(self, room: Room, time: float) -> float: """ Returns the rate at which air is being exchanged in the given room @@ -573,7 +570,24 @@ class ExposureModel: # Probability of infection. return (1 - np.exp(-inf_aero)) * 100 - def reproduction_rate(self): + def expected_new_cases(self): prob = self.infection_probability() exposed_occupants = self.exposed.number return prob * exposed_occupants / 100 + + def reproduction_number(self): + """ + The reproduction number can be thought of as the expected number of + cases directly generated by one infected case in a population. + + """ + if self.concentration_model.infected.number == 1: + return self.expected_new_cases() + + # Create an equivalent exposure model but with precisely + # one infected case. + single_exposure_model = nested_replace( + self, {'concentration_model.infected.number': 1} + ) + + return single_exposure_model.expected_new_cases() diff --git a/cara/tests/conftest.py b/cara/tests/conftest.py new file mode 100644 index 00000000..e50cdea8 --- /dev/null +++ b/cara/tests/conftest.py @@ -0,0 +1,38 @@ +from cara import models + +import pytest + + +@pytest.fixture +def baseline_model(): + 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'], + ), + ) + 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, + ) + ) diff --git a/cara/tests/test_dataclass_utils.py b/cara/tests/test_dataclass_utils.py new file mode 100644 index 00000000..ec0da382 --- /dev/null +++ b/cara/tests/test_dataclass_utils.py @@ -0,0 +1,27 @@ +import dataclasses + +from cara.dataclass_utils import nested_replace + + +@dataclasses.dataclass(frozen=True) +class Four: + four: float + + +@dataclasses.dataclass(frozen=True) +class Two: + three: int + four: Four + + +@dataclasses.dataclass(frozen=True) +class One: + one: int + two: Two + + + +def test_nested_replace(): + inst = One(1, two=Two(3, Four(4))) + new_inst = nested_replace(inst, {'two.four': Four(5)}) + assert new_inst == One(1, two=Two(3, Four(5))) diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index d621f31d..7118897b 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -24,41 +24,6 @@ def test_no_mask_emission_rate(baseline_model): ) -@pytest.fixture -def baseline_model(): - 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'], - ), - ) - 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( diff --git a/cara/tests/test_model.py b/cara/tests/test_model.py new file mode 100644 index 00000000..add6a899 --- /dev/null +++ b/cara/tests/test_model.py @@ -0,0 +1,12 @@ +import cara.models +from cara.dataclass_utils import nested_replace + + +def test_exposure_r0(baseline_exposure_model): + baseline_n3 = nested_replace( + baseline_exposure_model, {'concentration_model.infected.number': 3} + ) + # The number of new cases should be greater if there are more infecteds, but + # the reproduction number should be the same (it is a measure of one infected case). + assert baseline_n3.expected_new_cases() > baseline_exposure_model.expected_new_cases() + assert baseline_n3.reproduction_number() == baseline_exposure_model.reproduction_number()