diff --git a/cara/apps.py b/cara/apps.py index 09fd55f4..c6ee010e 100644 --- a/cara/apps.py +++ b/cara/apps.py @@ -49,13 +49,13 @@ class ConcentrationFigure: ax.set_xlabel('Time (hours)') ax.set_ylabel('Concentration ($q/m^3$)') ax.set_title('Concentration of infectious quanta aerosols') - ax.set_ymargin(0.2) - # ax.set_ylim(bottom=0) else: self.ax.ignore_existing_data_limits = True self.line.set_data(ts, concentration) - self.ax.relim() - self.ax.autoscale_view() + # Update the top limit based on the concentration if it exceeds 5 + # (rare but possible). + top = max([3, max(concentration)]) + self.ax.set_ylim(bottom=1e-4, top=top) self.figure.canvas.draw() @@ -160,14 +160,14 @@ class WidgetView: return widget def _build_window(self, node): - period = widgets.IntSlider(value=node.period, min=0, max=240) - interval = widgets.IntSlider(value=node.duration, min=0, max=240) + period = widgets.IntSlider(value=node.active.period, min=0, max=240) + interval = widgets.IntSlider(value=node.active.duration, min=0, max=240) def on_period_change(change): - node.period = change['new'] + node.active.period = change['new'] def on_interval_change(change): - node.duration = change['new'] + node.active.duration = change['new'] # TODO: Link the state back to the widget, not just the other way around. period.observe(on_period_change, names=['value']) @@ -181,14 +181,14 @@ class WidgetView: ) def _build_hepa(self, node): - period = widgets.IntSlider(value=node.period, min=0, max=240) - interval = widgets.IntSlider(value=node.duration, min=0, max=240) + period = widgets.IntSlider(value=node.active.period, min=0, max=240, step=5) + interval = widgets.IntSlider(value=node.active.duration, min=0, max=240, step=5) def on_period_change(change): - node.period = change['new'] + node.active.period = change['new'] def on_interval_change(change): - node.duration = change['new'] + node.active.duration = change['new'] # TODO: Link the state back to the widget, not just the other way around. period.observe(on_period_change, names=['value']) @@ -286,13 +286,14 @@ class WidgetView: baseline_model = models.Model( room=models.Room(volume=75), - ventilation=models.PeriodicWindow( - period=120, duration=120, inside_temp=293, outside_temp=283, cd_b=0.6, + ventilation=models.WindowOpening( + active=models.PeriodicInterval(period=120, duration=120), + inside_temp=293, outside_temp=283, cd_b=0.6, window_height=1.6, opening_length=0.6, ), infected=models.InfectedPerson( virus=models.Virus.types['SARS_CoV_2'], - present_times=((0, 4), (5, 8)), + 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'], @@ -313,14 +314,14 @@ class CARAStateBuilder(state.StateBuilder): def build_type_Ventilation(self, _: dataclasses.Field): s = state.DataclassStateNamed( states={ - 'Natural': self.build_generic(models.PeriodicWindow), - 'HEPA': self.build_generic(models.PeriodicHEPA), + 'Natural': self.build_generic(models.WindowOpening), + 'HEPA': self.build_generic(models.HEPAFilter), }, state_builder=self, ) # Initialise the HEPA state s._states['HEPA'].dcs_update_from( - models.PeriodicHEPA(120, 120, 500.) + models.HEPAFilter(models.PeriodicInterval(120, 120), 500.) ) return s diff --git a/cara/models.py b/cara/models.py index 998d025c..b74cdcf2 100644 --- a/cara/models.py +++ b/cara/models.py @@ -13,11 +13,79 @@ class Room: volume: int +@dataclass(frozen=True) +class Interval: + """ + Represents a collection of times in which a "thing" happens. + + The "thing" may be when an action is taken, such as opening a window, or + entering a room. + + Note that all intervals are open at the start, and closed at the end. So a + simple start, stop interval follows:: + + start < t <= end + + """ + def boundaries(self) -> typing.Tuple[typing.Tuple[float, float], ...]: + return () + + def transition_times(self) -> typing.Set[float]: + transitions = set() + for start, end in self.boundaries(): + transitions.update([start, end]) + return transitions + + def triggered(self, time: float) -> bool: + """Whether the given time falls inside this interval.""" + for start, end in self.boundaries(): + if start < time <= end: + return True + return False + + +@dataclass(frozen=True) +class SpecificInterval(Interval): + #: A sequence of times (start, stop), in hours, that the infected person + #: is present. The flattened list of times must be strictly monotonically + #: increasing. + present_times: typing.Tuple[typing.Tuple[float, float], ...] + + def boundaries(self): + return self.present_times + + +@dataclass(frozen=True) +class PeriodicInterval(Interval): + #: How often does the interval occur (minutes). + period: int + + #: How long does the interval occur for (minutes). + #: A value greater than :data:`period` signifies the event is permanently + #: occurring, a value of 0 signifies that the event never happens. + duration: int + + def boundaries(self) -> typing.Tuple[typing.Tuple[float, float], ...]: + result = [] + for i in np.arange(0, 24, self.period / 60): + result.append((i, i+self.duration/60)) + return tuple(result) + + @dataclass(frozen=True) class Ventilation: """ - An abstract class for ventilation schemes + Represents a mechanism by which air can be exchanged (replaced/filtered) + in a time dependent manner. + + The nature of the various air exchange schemes means that it is expected + for subclasses of Ventilation to exist. Known subclasses include + WindowOpening for window based air exchange, and HEPAFilter, for + mechanical air exchange through a filter. + """ + #: The times at which the air exchange is taking place. + active: Interval @abstractmethod def air_exchange(self, room: Room, time: float) -> float: @@ -25,16 +93,18 @@ class Ventilation: Returns the rate at which air is being exchanged in the given room per cubic meter at a given time (in hours). + Note that whilst the time is known inside this function, it may not + be used to vary the result unless the specific time used is declared + as part of a state change in the interval (e.g. when air_exchange == 0). + """ - pass + return 0. @dataclass(frozen=True) -class PeriodicWindow(Ventilation): - - # The window is opened for minutes every minutes - period: int #: How often the window is opened (minutes) - duration: int #: How long the window remains opened for (minutes) +class WindowOpening(Ventilation): + #: The interval in which the window is open. + active: Interval inside_temp: float #: The temperature inside the room (Kelvin) outside_temp: float #: The temperature outside of the window (Kelvin) @@ -46,37 +116,31 @@ class PeriodicWindow(Ventilation): cd_b: float = 0.6 #: Discharge coefficient: what portion effective area is used to exchange air (0 <= cd_b <= 1) def air_exchange(self, room: Room, time: float) -> float: - period = self.period / 60. - duration = self.duration / 60. - # Returns the rate at which air is being exchanged in the given room per cubic meter at a given time - # If the window is closed, no air is being exchanged - if (time % period) < (period - duration): - return 0 + # If the window is shut, no air is being exchanged. + if not self.active.triggered(time): + return 0. - root = np.sqrt(9.81 * self.window_height * (abs(self.inside_temp - self.outside_temp)) / self.outside_temp) + # Reminder, no dependence on time in the resulting calculation. + + temp_delta = abs(self.inside_temp - self.outside_temp) / self.outside_temp + root = np.sqrt(9.81 * self.window_height * temp_delta) return (3600 / (3 * room.volume)) * self.cd_b * self.window_height * self.opening_length * root @dataclass(frozen=True) -class PeriodicHEPA(Ventilation): +class HEPAFilter(Ventilation): + #: The interval in which the HEPA filter is operating. + active: Interval - # The HEPA is switched on for minutes every minutes - period: int #: How often the HEPA is switched on (minutes) - duration: int #: How long the HEPA remains switched on for (minutes) - - q_air_mech: float #: The rate at which the HEPA exchanges air (when switched on) + #: The rate at which the HEPA exchanges air (when switched on) + q_air_mech: float def air_exchange(self, room: Room, time: float) -> float: - # Returns the rate at which air is being exchanged in the given room per cubic meter at a given time - - period = self.period / 60. - duration = self.duration / 60. - - # If the HEPA is off, no air is being exchanged - if (time % period) < (period - duration): - return 0 - + # If the HEPA is off, no air is being exchanged. + if not self.active.triggered(time): + return 0. + # Reminder, no dependence on time in the resulting calculation. return self.q_air_mech / room.volume @@ -197,30 +261,14 @@ Activity.types = { @dataclass(frozen=True) class InfectedPerson: virus: Virus - #: A sequence of times (start, stop), in hours, that the infected person - #: is present. The flattened list of times must be strictly monotonically - #: increasing. - present_times: typing.Tuple[typing.Tuple[float, float], ...] + #: The times in which the person is in the room. + presence: Interval mask: Mask activity: Activity expiration: Expiration def person_present(self, time): - for start, end in self.present_times: - if start <= time <= end: - return True - return False - - def start_end_of_presence(self, time) -> typing.Tuple[float, float]: - """ - Find the most recent start (and associated) end-time (even if the - given time is after the end-point) given a time. - - """ - for start, end in self.present_times[::-1]: - if time > start: - return start, end - return start, end + return self.presence.triggered(time) @functools.lru_cache() def emission_rate(self, time) -> float: @@ -266,34 +314,44 @@ class Model: return k + self.virus.decay_constant + self.ventilation.air_exchange(self.room, time) - def collect_time_state_changes(self): + @functools.lru_cache() + def state_change_times(self): """ All time dependent entities on this model must provide information about the times at which their state changes. """ + state_change_times = set() + state_change_times.update(self.infected.presence.transition_times()) + state_change_times.update(self.ventilation.active.transition_times()) + return sorted(state_change_times) + + def last_state_change(self, time: float): + """ + Find the most recent state change. + + """ + for change_time in self.state_change_times()[::-1]: + if change_time < time: + return change_time + return 0 @functools.lru_cache() def concentration(self, time: float) -> float: - t = time + if time == 0: + return 0.0 IVRR = self.infectious_virus_removal_rate(time) V = self.room.volume Ni = self.infected_occupants ER = self.infected.emission_rate(time) - t0, t1 = self.infected.start_end_of_presence(time) - if t == 0: - return 0.0 - elif t0 < t <= t1: - # Concentration while infected present. - init_concentration = self.concentration(t0) - fac = np.exp(IVRR * (t0 - t)) - return ((ER * Ni) / (IVRR * V)) * (1 - fac) + init_concentration * fac - else: - # Concentration while infected not present. - end_concentration = self.concentration(t1) - fac = np.exp(IVRR * (t1 - t)) - return end_concentration * fac + 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) + return concentration_limit * (1 - fac) + concentration_at_last_state_change * fac def infection_probability(self): # Infection probability @@ -306,7 +364,7 @@ class Model: return np.trapz([fn(v) for v in values], values) # TODO: Have this for exposed not infected. - for start, stop in self.infected.present_times: + for start, stop in self.infected.presence.boundaries(): exposure += (integrate(self.concentration, start, stop)) inf_aero = ( @@ -317,4 +375,3 @@ class Model: # Probability of infection. return (1 - np.exp(-inf_aero)) * 100 - diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py index 9dc36f38..5a11e4ae 100644 --- a/cara/tests/test_known_quantities.py +++ b/cara/tests/test_known_quantities.py @@ -1,3 +1,4 @@ +import numpy as np import numpy.testing as npt import pytest @@ -17,7 +18,7 @@ def test_no_mask_emission_rate(baseline_model): rate = 167.74011998223307 npt.assert_allclose( [baseline_model.infected.emission_rate(t) for t in [0, 1, 4, 4.5, 5, 8, 9]], - [rate, rate, rate, 0, rate, rate, 0], + [0, rate, rate, 0, 0, rate, 0], rtol=1e-5 ) @@ -26,11 +27,14 @@ def test_no_mask_emission_rate(baseline_model): def baseline_model(): model = models.Model( room=models.Room(volume=75), - ventilation=models.PeriodicWindow(period=120, duration=120, inside_temp=293, outside_temp=283, cd_b=0.6, - window_height=1.6, opening_length=0.6), + ventilation=models.WindowOpening( + active=models.PeriodicInterval(period=120, duration=120), + inside_temp=293, outside_temp=283, cd_b=0.6, + window_height=1.6, opening_length=0.6, + ), infected=models.InfectedPerson( virus=models.Virus.types['SARS_CoV_2'], - present_times=((0, 4), (5, 8)), + 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'], @@ -44,8 +48,11 @@ def baseline_model(): @pytest.fixture def baseline_periodic_window(): - return models.PeriodicWindow(period=120, duration=15, inside_temp=293, outside_temp=283, cd_b=0.6, - window_height=1.6, opening_length=0.6) + return models.WindowOpening( + active=models.PeriodicInterval(period=120, duration=15), + inside_temp=293, outside_temp=283, cd_b=0.6, + window_height=1.6, opening_length=0.6, + ) @pytest.fixture @@ -55,7 +62,10 @@ def baseline_room(): @pytest.fixture def baseline_periodic_hepa(): - return models.PeriodicHEPA(period=120, duration=15, q_air_mech=514.74) + return models.HEPAFilter( + active=models.PeriodicInterval(period=120, duration=15), + q_air_mech=514.74, + ) def test_concentrations(baseline_model): @@ -68,6 +78,46 @@ def test_concentrations(baseline_model): ) +def test_smooth_concentrations(baseline_model): + # We don't care about the actual concentrations in this test, but rather + # that the curve itself is smooth. + dx = 0.1 + dy_limit = dx * 2 # Anything more than this is a bit steep. + ts = np.arange(0, 10, dx) + concentrations = [baseline_model.concentration(t) for t in ts] + assert np.abs(np.diff(concentrations)).max() < dy_limit + + +def build_model(interval_duration): + model = models.Model( + room=models.Room(volume=75), + ventilation=models.HEPAFilter( + active=models.PeriodicInterval(period=120, duration=interval_duration), + q_air_mech=500., + ), + infected=models.InfectedPerson( + 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 + + +def test_concentrations_startup(baseline_model): + # The concentrations should be the same until the beginning of the + # first time that the ventilation is disabled. + m1 = build_model(interval_duration=120) + m2 = build_model(interval_duration=65) + + assert m1.concentration(1.) == m2.concentration(1.) + + def test_r0(baseline_model): p = baseline_model.infection_probability() npt.assert_allclose(p, 93.196908) @@ -78,7 +128,7 @@ def test_periodic_window(baseline_periodic_window, baseline_room): ts = [0, 14/60, 15/60, 16/60, 119/60, 120/60, 121/60, 239/60, 240/60] aes = [baseline_periodic_window.air_exchange(baseline_room, t) for t in ts] rate = 6.86347 - answers = [0, 0, 0, 0, rate, 0, 0, rate, 0] + answers = [0, rate, rate, 0, 0, 0, rate, 0, 0] npt.assert_allclose(aes, answers, rtol=1e-5) @@ -87,5 +137,5 @@ def test_periodic_hepa(baseline_periodic_hepa, baseline_room): ts = [0, 14 / 60, 15 / 60, 16 / 60, 119 / 60, 120 / 60, 121 / 60, 239 / 60, 240 / 60] rate = 514.74 / 75 aes = [baseline_periodic_hepa.air_exchange(baseline_room, t) for t in ts] - answers = [0, 0, 0, 0, rate, 0, 0, rate, 0] + answers = [0, rate, rate, 0, 0, 0, rate, 0, 0] npt.assert_allclose(aes, answers, rtol=1e-5)