From a70f899b5238f09d529a87edcc7dc24c0a1be7d1 Mon Sep 17 00:00:00 2001
From: Phil Elson <philip.elson@cern.ch>
Date: Tue, 27 Oct 2020 06:27:38 +0100
Subject: [PATCH 1/4] Extract the concept of an interval for the Ventilation
 schemes, and correctly handle the time dependence of Ventilation in the
 concentration calculations.

---
 cara/apps.py                        |  27 +++---
 cara/models.py                      | 144 ++++++++++++++++++++--------
 cara/tests/test_known_quantities.py |  21 ++--
 3 files changed, 134 insertions(+), 58 deletions(-)

diff --git a/cara/apps.py b/cara/apps.py
index 09fd55f4..2fad786c 100644
--- a/cara/apps.py
+++ b/cara/apps.py
@@ -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,8 +286,9 @@ 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(
@@ -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..ce0c64d9 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -13,11 +13,68 @@ 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.
+
+    """
+    def triggered(self, time: float) -> bool:
+        """Whether the given time falls inside this interval."""
+        return False
+
+    def boundaries(self) -> typing.Set[float]:
+        """Returns the edges of this interval."""
+        return set()
+
+
+@dataclass(frozen=True)
+class PeriodicInterval(Interval):
+    #: How often does the interval occur.
+    period: int
+
+    #: How long does the interval occur for.
+    #: 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 triggered(self, time: float) -> bool:
+        period = self.period / 60.
+        duration = self.duration / 60.
+        return (time % period) >= (period - duration)
+
+    def boundaries(self) -> typing.Set[float]:
+        state_changes = set()
+        period_h = self.period / 60
+        duration_h = self.duration / 60
+        # The current implementation starts at the end of the first period, not at 0.
+        start = period_h
+        # Take as many steps as we need to get to a full day.
+        for i in np.arange(start, 24, period_h):
+            state_changes.add(i)
+            if duration_h < period_h:
+                state_changes.add(i - duration_h)
+
+        return state_changes
+
+
 @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:
@@ -26,15 +83,20 @@ class Ventilation:
         cubic meter at a given time (in hours).
 
         """
-        pass
+        return 0.
+
+    def times_of_state_change(self) -> typing.Set[float]:
+        """
+        Returns the times at which a change in ventilation occurs.
+
+        """
+        return self.active.boundaries()
 
 
 @dataclass(frozen=True)
-class PeriodicWindow(Ventilation):
-
-    # The window is opened for <duration> minutes every <period> 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,36 +108,27 @@ 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)
+        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):
-
-    # The HEPA is switched on for <duration> minutes every <period> minutes
-    period: int   #: How often the HEPA is switched on (minutes)
-    duration: int   #: How long the HEPA remains switched on for (minutes)
+class HEPAFilter(Ventilation):
+    #: The interval in which the HEPA filter is operating.
+    active: Interval
 
     q_air_mech: float   #: The rate at which the HEPA exchanges air (when switched on)
 
     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.
 
         return self.q_air_mech / room.volume
 
@@ -266,12 +319,29 @@ class Model:
 
         return k + self.virus.decay_constant + self.ventilation.air_exchange(self.room, time)
 
+    @functools.lru_cache()
     def collect_time_state_changes(self):
         """
         All time dependent entities on this model must provide information about
         the times at which their state changes.
 
         """
+        state_change_times = set()
+        for start, end in self.infected.present_times:
+            state_change_times.add(start)
+            state_change_times.add(end)
+        state_change_times.update(self.ventilation.times_of_state_change())
+        return sorted(state_change_times)
+
+    def last_state_change(self, time: float):
+        """
+        Find the most recent state change.
+
+        """
+        for change_time in self.collect_time_state_changes()[::-1]:
+            if time >= change_time:
+                return change_time
+        return 0
 
     @functools.lru_cache()
     def concentration(self, time: float) -> float:
@@ -280,20 +350,17 @@ class Model:
         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
+        t_last_state_change = self.last_state_change(time)
+        if t != t_last_state_change:
+            concentration_at_last_state_change = self.concentration(t_last_state_change)
         else:
-            # Concentration while infected not present.
-            end_concentration = self.concentration(t1)
-            fac = np.exp(IVRR * (t1 - t))
-            return end_concentration * fac
+            concentration_at_last_state_change = self.concentration(t - 0.0001)
+
+        fac = np.exp(IVRR * (t_last_state_change - t))
+        concentration_limit = (ER * Ni) / (IVRR * V)
+        return concentration_limit * (1 - fac) + concentration_at_last_state_change * fac
 
     def infection_probability(self):
         # Infection probability
@@ -317,4 +384,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..b6e4be51 100644
--- a/cara/tests/test_known_quantities.py
+++ b/cara/tests/test_known_quantities.py
@@ -26,8 +26,11 @@ 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)),
@@ -44,8 +47,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 +61,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):
@@ -63,7 +72,7 @@ def test_concentrations(baseline_model):
     concentrations = [baseline_model.concentration(t) for t in ts]
     npt.assert_allclose(
         concentrations,
-        [0.000000e+00, 2.891970e-01, 1.266287e-04, 2.891969e-01, 5.544607e-08],
+        [0.000000e+00, 2.891970e-01, 1.267267e-04, 2.891969e-01, 5.548897e-08],
         rtol=1e-5
     )
 

From 2a20bf834f31657e05a716bb5a3409c30faa5239 Mon Sep 17 00:00:00 2001
From: Phil Elson <philip.elson@cern.ch>
Date: Tue, 27 Oct 2020 14:34:45 +0100
Subject: [PATCH 2/4] Improve the concentration calculation scheme by more
 strictly defining the time interval concept to be closed start and open end

---
 cara/apps.py                        |  9 +++-
 cara/models.py                      | 81 ++++++++++++++++++-----------
 cara/tests/test_known_quantities.py | 49 +++++++++++++++--
 3 files changed, 103 insertions(+), 36 deletions(-)

diff --git a/cara/apps.py b/cara/apps.py
index 2fad786c..3caa3fb1 100644
--- a/cara/apps.py
+++ b/cara/apps.py
@@ -34,8 +34,8 @@ class ConcentrationFigure:
         self.line = None
 
     def update(self, model: models.Model):
-        resolution = 600
-        ts = np.linspace(0, 10, resolution)
+        resolution = 500
+        ts = np.linspace(0, 9.09, resolution)
         concentration = [model.concentration(t) for t in ts]
         if self.line is None:
             [self.line] = self.ax.plot(ts, concentration)
@@ -56,6 +56,11 @@ class ConcentrationFigure:
             self.line.set_data(ts, concentration)
             self.ax.relim()
             self.ax.autoscale_view()
+        # self.ax.set_yscale('log')
+        # if max(concentration) > 1:
+        self.ax.set_ylim(bottom=1e-4, top=5)
+        # else:
+        #     self.ax.set_ylim(bottom=0, top=1)
         self.figure.canvas.draw()
 
 
diff --git a/cara/models.py b/cara/models.py
index ce0c64d9..643d33ea 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -33,31 +33,54 @@ class Interval:
 
 @dataclass(frozen=True)
 class PeriodicInterval(Interval):
-    #: How often does the interval occur.
+    #: How often does the interval occur (minutes).
     period: int
 
-    #: How long does the interval occur for.
+    #: 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 states(self):
+        result = []
+        for i in np.arange(0, 24, self.period / 60):
+            result.append((i, i+self.duration/60))
+        return tuple(result)
+
     def triggered(self, time: float) -> bool:
+        for start, end in self.states():
+            if start < time <= end:
+                return True
+        return False
+
+        if self.duration >= self.period:
+            return True
+        if self.duration == 0:
+            return False
+
         period = self.period / 60.
         duration = self.duration / 60.
-        return (time % period) >= (period - duration)
+        t = time % period
+        return 0 <= t < duration
 
     def boundaries(self) -> typing.Set[float]:
         state_changes = set()
+        for start, end in self.states():
+            state_changes.add(start)
+            state_changes.add(end)
+        return state_changes
+
+        if self.duration >= self.period:
+            return set()
+        if self.duration == 0:
+            return set()
+        state_changes = set()
         period_h = self.period / 60
         duration_h = self.duration / 60
-        # The current implementation starts at the end of the first period, not at 0.
-        start = period_h
         # Take as many steps as we need to get to a full day.
-        for i in np.arange(start, 24, period_h):
+        for i in np.arange(0, 24, period_h):
             state_changes.add(i)
-            if duration_h < period_h:
-                state_changes.add(i - duration_h)
-
+            state_changes.add(i + duration_h)
         return state_changes
 
 
@@ -260,20 +283,20 @@ class InfectedPerson:
 
     def person_present(self, time):
         for start, end in self.present_times:
-            if start <= time <= end:
+            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
+    # 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
 
     @functools.lru_cache()
     def emission_rate(self, time) -> float:
@@ -339,26 +362,24 @@ class Model:
 
         """
         for change_time in self.collect_time_state_changes()[::-1]:
-            if time >= change_time:
+            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)
-        if t == 0:
-            return 0.0
-        t_last_state_change = self.last_state_change(time)
-        if t != t_last_state_change:
-            concentration_at_last_state_change = self.concentration(t_last_state_change)
-        else:
-            concentration_at_last_state_change = self.concentration(t - 0.0001)
 
-        fac = np.exp(IVRR * (t_last_state_change - t))
+        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
 
diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py
index b6e4be51..82393d4c 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
     )
 
@@ -72,11 +73,51 @@ def test_concentrations(baseline_model):
     concentrations = [baseline_model.concentration(t) for t in ts]
     npt.assert_allclose(
         concentrations,
-        [0.000000e+00, 2.891970e-01, 1.267267e-04, 2.891969e-01, 5.548897e-08],
+        [0.000000e+00, 2.891970e-01, 1.266287e-04, 2.891969e-01, 5.544607e-08],
         rtol=1e-5
     )
 
 
+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'],
+            present_times=((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)
@@ -87,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)
 
 
@@ -96,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)

From 332f2414ad16591b97b692d111e7ca82756d040e Mon Sep 17 00:00:00 2001
From: Phil Elson <philip.elson@cern.ch>
Date: Tue, 27 Oct 2020 14:47:45 +0100
Subject: [PATCH 3/4] Abstract the presence time as an interval.

---
 cara/apps.py                        |  6 +--
 cara/models.py                      | 80 ++++++++++++-----------------
 cara/tests/test_known_quantities.py |  4 +-
 3 files changed, 39 insertions(+), 51 deletions(-)

diff --git a/cara/apps.py b/cara/apps.py
index 3caa3fb1..dbca65fb 100644
--- a/cara/apps.py
+++ b/cara/apps.py
@@ -34,8 +34,8 @@ class ConcentrationFigure:
         self.line = None
 
     def update(self, model: models.Model):
-        resolution = 500
-        ts = np.linspace(0, 9.09, resolution)
+        resolution = 600
+        ts = np.linspace(0, 10, resolution)
         concentration = [model.concentration(t) for t in ts]
         if self.line is None:
             [self.line] = self.ax.plot(ts, concentration)
@@ -298,7 +298,7 @@ baseline_model = models.Model(
     ),
     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'],
diff --git a/cara/models.py b/cara/models.py
index 643d33ea..07b7f382 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -21,6 +21,11 @@ class Interval:
     The "thing" may be when an action is taken, such as opening a window, or
     entering a room.
 
+    Note that all intervals are closed at the start, and open at the end. So a
+    simple start, stop interval follows::
+
+        start < t <= end
+
     """
     def triggered(self, time: float) -> bool:
         """Whether the given time falls inside this interval."""
@@ -31,6 +36,27 @@ class Interval:
         return set()
 
 
+@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 triggered(self, time: float) -> bool:
+        for start, end in self.present_times:
+            if start < time <= end:
+                return True
+        return False
+
+    def boundaries(self) -> typing.Set[float]:
+        state_changes = set()
+        for start, end in self.present_times:
+            state_changes.add(start)
+            state_changes.add(end)
+        return state_changes
+
+
 @dataclass(frozen=True)
 class PeriodicInterval(Interval):
     #: How often does the interval occur (minutes).
@@ -53,16 +79,6 @@ class PeriodicInterval(Interval):
                 return True
         return False
 
-        if self.duration >= self.period:
-            return True
-        if self.duration == 0:
-            return False
-
-        period = self.period / 60.
-        duration = self.duration / 60.
-        t = time % period
-        return 0 <= t < duration
-
     def boundaries(self) -> typing.Set[float]:
         state_changes = set()
         for start, end in self.states():
@@ -70,19 +86,6 @@ class PeriodicInterval(Interval):
             state_changes.add(end)
         return state_changes
 
-        if self.duration >= self.period:
-            return set()
-        if self.duration == 0:
-            return set()
-        state_changes = set()
-        period_h = self.period / 60
-        duration_h = self.duration / 60
-        # Take as many steps as we need to get to a full day.
-        for i in np.arange(0, 24, period_h):
-            state_changes.add(i)
-            state_changes.add(i + duration_h)
-        return state_changes
-
 
 @dataclass(frozen=True)
 class Ventilation:
@@ -273,30 +276,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:
@@ -350,9 +337,10 @@ class Model:
 
         """
         state_change_times = set()
-        for start, end in self.infected.present_times:
-            state_change_times.add(start)
-            state_change_times.add(end)
+        # for start, end in self.infected.present_times:
+        #     state_change_times.add(start)
+        #     state_change_times.add(end)
+        state_change_times.update(self.infected.presence.boundaries())
         state_change_times.update(self.ventilation.times_of_state_change())
         return sorted(state_change_times)
 
@@ -394,7 +382,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.present_times:
             exposure += (integrate(self.concentration, start, stop))
 
         inf_aero = (
diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py
index 82393d4c..5a11e4ae 100644
--- a/cara/tests/test_known_quantities.py
+++ b/cara/tests/test_known_quantities.py
@@ -34,7 +34,7 @@ def baseline_model():
         ),
         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'],
@@ -97,7 +97,7 @@ def build_model(interval_duration):
         ),
         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'],

From b8b322849b6fdcf523a99cb7bb66573f828505c2 Mon Sep 17 00:00:00 2001
From: Phil Elson <philip.elson@cern.ch>
Date: Tue, 27 Oct 2020 15:06:28 +0100
Subject: [PATCH 4/4] Code review.

---
 cara/apps.py   | 13 +++------
 cara/models.py | 78 +++++++++++++++++++-------------------------------
 2 files changed, 34 insertions(+), 57 deletions(-)

diff --git a/cara/apps.py b/cara/apps.py
index dbca65fb..c6ee010e 100644
--- a/cara/apps.py
+++ b/cara/apps.py
@@ -49,18 +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()
-        # self.ax.set_yscale('log')
-        # if max(concentration) > 1:
-        self.ax.set_ylim(bottom=1e-4, top=5)
-        # else:
-        #     self.ax.set_ylim(bottom=0, top=1)
+        # 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()
 
 
diff --git a/cara/models.py b/cara/models.py
index 07b7f382..b74cdcf2 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -21,20 +21,28 @@ class Interval:
     The "thing" may be when an action is taken, such as opening a window, or
     entering a room.
 
-    Note that all intervals are closed at the start, and open at the end. So a
+    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
 
-    def boundaries(self) -> typing.Set[float]:
-        """Returns the edges of this interval."""
-        return set()
-
 
 @dataclass(frozen=True)
 class SpecificInterval(Interval):
@@ -43,18 +51,8 @@ class SpecificInterval(Interval):
     #: increasing.
     present_times: typing.Tuple[typing.Tuple[float, float], ...]
 
-    def triggered(self, time: float) -> bool:
-        for start, end in self.present_times:
-            if start < time <= end:
-                return True
-        return False
-
-    def boundaries(self) -> typing.Set[float]:
-        state_changes = set()
-        for start, end in self.present_times:
-            state_changes.add(start)
-            state_changes.add(end)
-        return state_changes
+    def boundaries(self):
+        return self.present_times
 
 
 @dataclass(frozen=True)
@@ -67,25 +65,12 @@ class PeriodicInterval(Interval):
     #: occurring, a value of 0 signifies that the event never happens.
     duration: int
 
-    def states(self):
+    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)
 
-    def triggered(self, time: float) -> bool:
-        for start, end in self.states():
-            if start < time <= end:
-                return True
-        return False
-
-    def boundaries(self) -> typing.Set[float]:
-        state_changes = set()
-        for start, end in self.states():
-            state_changes.add(start)
-            state_changes.add(end)
-        return state_changes
-
 
 @dataclass(frozen=True)
 class Ventilation:
@@ -108,16 +93,13 @@ 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).
+
         """
         return 0.
 
-    def times_of_state_change(self) -> typing.Set[float]:
-        """
-        Returns the times at which a change in ventilation occurs.
-
-        """
-        return self.active.boundaries()
-
 
 @dataclass(frozen=True)
 class WindowOpening(Ventilation):
@@ -138,6 +120,8 @@ class WindowOpening(Ventilation):
         if not self.active.triggered(time):
             return 0.
 
+        # 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)
 
@@ -149,13 +133,14 @@ class HEPAFilter(Ventilation):
     #: The interval in which the HEPA filter is operating.
     active: Interval
 
-    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:
         # 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
 
 
@@ -330,18 +315,15 @@ class Model:
         return k + self.virus.decay_constant + self.ventilation.air_exchange(self.room, time)
 
     @functools.lru_cache()
-    def collect_time_state_changes(self):
+    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()
-        # for start, end in self.infected.present_times:
-        #     state_change_times.add(start)
-        #     state_change_times.add(end)
-        state_change_times.update(self.infected.presence.boundaries())
-        state_change_times.update(self.ventilation.times_of_state_change())
+        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):
@@ -349,7 +331,7 @@ class Model:
         Find the most recent state change.
 
         """
-        for change_time in self.collect_time_state_changes()[::-1]:
+        for change_time in self.state_change_times()[::-1]:
             if change_time < time:
                 return change_time
         return 0
@@ -382,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.presence.present_times:
+        for start, stop in self.infected.presence.boundaries():
             exposure += (integrate(self.concentration, start, stop))
 
         inf_aero = (