Recording the serial interval

pyEpiabm/pyEpiabm/core/person.py

@@ -54,6 +54,8 @@ def __init__(self, microcell, age_group=None):
         self.secondary_infections_counts = []
         self.time_of_recovery = None
         self.num_times_infected = 0
+        self.exposure_period = None
+        self.serial_interval_dict = {}
         self.care_home_resident = False
         self.key_worker = False
         self.date_positive = None
@@ -283,10 +285,56 @@ def increment_secondary_infections(self):
         """Increments the number of secondary infections the given person has
         for this specific infection period (i.e. if the given person has been
         infected multiple times, then we only increment the current secondary
-        infection count)
+        infection count).
         """
         try:
             self.secondary_infections_counts[-1] += 1
         except IndexError:
             raise RuntimeError("Cannot call increment_secondary_infections "
                                "while secondary_infections_counts is empty")
+
+    def set_exposure_period(self, exposure_period: float):
+        """Sets the exposure period (we define here as the time between a
+        primary case infection and a secondary case exposure, with the current
+        `Person` being the secondary case). We store this to be added to the
+        latent period of the infection to give a serial interval.
+
+        Parameters
+        ----------
+        exposure_period : float
+            The time between the infector's time of infection and the time
+            of exposure to the current Person
+        """
+        self.exposure_period = exposure_period
+
+    def store_serial_interval(self, latent_period: float):
+        """Adds this `latent_period` to the current `exposure_period` to give
+        a `serial_interval`, which will be stored in the
+        `serial_interval_dict`. The serial interval is the time between a
+        primary case infection and a secondary case infection. This method
+        is called immediately after a person becomes exposed.
+
+        Parameters
+        ----------
+        latent_period : float
+            The period between this `Person`'s time of exposure and this
+            `Person`'s time of infection
+        """
+        # This method has been called erroneously if the exposure period is
+        # None
+        if self.exposure_period is None:
+            raise RuntimeError("Cannot call store_serial_interval while the"
+                               " exposure_period is None")
+
+        serial_interval = self.exposure_period + latent_period
+        # The reference day is the day the primary case was first infected
+        # This is what we will store in the dictionary
+        reference_day = self.time_of_status_change - serial_interval
+        try:
+            (self.serial_interval_dict[reference_day]
+             .append(serial_interval))
+        except KeyError:
+            self.serial_interval_dict[reference_day] = [serial_interval]
+
+        # Reset the exposure period for the next infection
+        self.exposure_period = None

pyEpiabm/pyEpiabm/routine/simulation.py

@@ -65,6 +65,8 @@ def configure(self,
                a csv file containing infectiousness (viral load) values
             * `secondary_infections_output`: Boolean to determine whether we \
                need a csv file containing secondary infections and R_t values
+            * `serial_interval_output`: Boolean to determine whether we \
+               need a csv file containing serial interval data
             * `compress`: Boolean to determine whether we compress \
                the infection history csv files
 
@@ -87,12 +89,12 @@ def configure(self,
         inf_history_params : dict
             This is short for 'infection history file parameters' and we will
             use the abbreviation 'ih' to refer to infection history throughout
-            this class. If `status_output`, `infectiousness_output` and
-            `secondary_infections_output` are all False, then no infection
-            history csv files are produced (or if the dictionary is None).
-            These files contain the infection status, infectiousness and
-            secondary infection counts of each person every time step
-            respectively. The EpiOS tool
+            this class. If `status_output`, `infectiousness_output`,
+            `secondary_infections_output` and `serial_interval_output` are all
+            False, then no infection history csv files are produced (or if the
+            dictionary is None). These files contain the infection status,
+            infectiousness and secondary infection counts of each person every
+            time step respectively. The EpiOS tool
             (https://github.com/SABS-R3-Epidemiology/EpiOS) samples data from
             these files to mimic real life epidemic sampling techniques. These
             files can be compressed when 'compress' is True, reducing the size
@@ -154,9 +156,11 @@ def configure(self,
         self.status_output = False
         self.infectiousness_output = False
         self.secondary_infections_output = False
+        self.serial_interval_output = False
         self.ih_status_writer = None
         self.ih_infectiousness_writer = None
         self.secondary_infections_writer = None
+        self.serial_interval_writer = None
         self.compress = False
 
         if inf_history_params:
@@ -168,19 +172,28 @@ def configure(self,
                 .get("infectiousness_output")
             self.secondary_infections_output = inf_history_params\
                 .get("secondary_infections_output")
+            self.serial_interval_output = inf_history_params \
+                .get("serial_interval_output")
             self.compress = inf_history_params.get("compress", False)
             person_ids = []
             person_ids += [person.id for cell in population.cells for person
                            in cell.persons]
             self.ih_output_titles = ["time"] + person_ids
             self.Rt_output_titles = ["time"] + person_ids + ["R_t"]
+            ts = 1 / Parameters.instance().time_steps_per_day
+            times = np.arange(self.sim_params["simulation_start_time"],
+                              self.sim_params["simulation_end_time"] + ts,
+                              ts).tolist()
+            self.si_output_titles = times
             ih_folder = os.path.join(os.getcwd(),
                                      inf_history_params["output_dir"])
 
             if not (self.status_output or self.infectiousness_output
-                    or self.secondary_infections_output):
-                logging.warning("status_output, infectiousness_output and "
-                                + "secondary_infections_output are False. "
+                    or self.secondary_infections_output
+                    or self.serial_interval_output):
+                logging.warning("status_output, infectiousness_output, "
+                                + "secondary_infections_output and "
+                                + "serial_interval_output are False. "
                                 + "No infection history csvs will be created.")
 
             if self.status_output:
@@ -219,6 +232,18 @@ def configure(self,
                     self.Rt_output_titles
                 )
 
+            if self.serial_interval_output:
+
+                file_name = "serial_intervals.csv"
+                logging.info(
+                    f"Set serial interval location to "
+                    f"{os.path.join(ih_folder, file_name)}")
+
+                self.serial_interval_writer = _CsvDictWriter(
+                    ih_folder, file_name,
+                    self.si_output_titles
+                )
+
     @log_exceptions()
     def run_sweeps(self):
         """Iteration step of the simulation. First the initialisation sweeps
@@ -264,6 +289,8 @@ def run_sweeps(self):
         logging.info(f"Final time {t} days reached")
         if self.secondary_infections_writer:
             self.write_to_Rt_file(times)
+        if self.serial_interval_writer:
+            self.write_to_serial_interval_file(times)
 
     def write_to_file(self, time):
         """Records the count number of a given list of infection statuses
@@ -404,6 +431,46 @@ def write_to_Rt_file(self, times: np.array):
             # Write each time step in dictionary form
             self.secondary_infections_writer.write(dict_row)
 
+    def write_to_serial_interval_file(self, times: np.array):
+        """Records the intervals between an infector and an infectee getting
+        infected to provide an overall serial interval for each time-step of
+        the epidemic. This can be used as a histogram of values for each
+        time step.
+
+        Parameters
+        ----------
+        times : np.array
+            An array of all time steps of the simulation
+        """
+        # Initialise the dataframe
+        all_times = np.hstack((np.array(self
+                                        .sim_params["simulation_start_time"]),
+                               times))
+        data_dict = {time: [] for time in all_times}
+        for cell in self.population.cells:
+            for person in cell.persons:
+                # For every time the person was infected, add their list of
+                # serial intervals to the timepoint at which their infector
+                # became infected
+                for t_inf, intervals in person.serial_interval_dict.items():
+                    data_dict[t_inf] += intervals
+
+        # Here we will fill out the rest of the dataframe with NaN values,
+        # as all lists will have different lengths
+        max_list_length = max([len(intervals)
+                               for intervals in data_dict.values()])
+        for t in data_dict.keys():
+            data_dict[t] += ([np.nan] * (max_list_length - len(data_dict[t])))
+
+        # Change to dataframe to get the data in a list of dicts format
+        df = pd.DataFrame(data_dict)
+
+        # The below is a list of dictionaries for each time step
+        list_of_dicts = df.to_dict(orient='records')
+        for dict_row in list_of_dicts:
+            # Write each time step in dictionary form
+            self.serial_interval_writer.write(dict_row)
+
     def add_writer(self, writer: AbstractReporter):
         self.writers.append(writer)
 

pyEpiabm/pyEpiabm/sweep/host_progression_sweep.py

@@ -324,6 +324,13 @@ def update_time_status_change(self, person: Person, time: float):
 
         person.time_of_status_change = time + transition_time
 
+        # Finally, if the person is Exposed, we can store their latency period
+        # as the transition_time. This can be used for calculating the serial
+        # interval
+        if person.infection_status == InfectionStatus.Exposed:
+            latent_period = transition_time
+            person.store_serial_interval(latent_period)
+
     def _updates_infectiousness(self, person: Person, time: float):
         """Updates infectiousness. Scales using the initial infectiousness
         if the person is in an infectious state. Updates the infectiousness to

pyEpiabm/pyEpiabm/sweep/household_sweep.py

@@ -56,3 +56,9 @@ def __call__(self, time: float):
                         # Increment the infector's
                         # secondary_infections_count
                         infector.increment_secondary_infections()
+                        # Set the time between infector's infection time and
+                        # the infectee's exposure time (current time) to be
+                        # the exposure period of the infectee
+                        inf_to_exposed = (time -
+                                          infector.infection_start_times[-1])
+                        infectee.set_exposure_period(inf_to_exposed)

pyEpiabm/pyEpiabm/sweep/place_sweep.py

@@ -57,6 +57,12 @@ def __call__(self, time: float):
                             # Increment the infector's
                             # secondary_infections_count
                             infector.increment_secondary_infections()
+                            # Set the time between infector's infection time
+                            # and the infectee's exposure time (current time)
+                            # to be the exposure period of the infectee
+                            inf_to_exposed = \
+                                (time - infector.infection_start_times[-1])
+                            infectee.set_exposure_period(inf_to_exposed)
 
                     # Otherwise number of infectees is binomially
                     # distributed. Not sure if covidsim considers only
@@ -97,3 +103,10 @@ def __call__(self, time: float):
                                 # Increment the infector's
                                 # secondary_infections_count
                                 infector.increment_secondary_infections()
+                                # Set the time between infector's infection
+                                # time and the infectee's exposure time
+                                # (current time) to be the exposure period of
+                                # the infectee
+                                inf_to_exposed = \
+                                    (time - infector.infection_start_times[-1])
+                                infectee.set_exposure_period(inf_to_exposed)

pyEpiabm/pyEpiabm/sweep/spatial_sweep.py

@@ -251,6 +251,12 @@ def do_infection_event(self, infector: Person, infectee: Person,
             infectee.microcell.cell.enqueue_person(infectee)
             # Increment the infector's secondary_infections_count
             infector.increment_secondary_infections()
+            # Set the time between infector's infection time and
+            # the infectee's exposure time (current time) to be
+            # the exposure period of the infectee
+            inf_to_exposed = (time -
+                              infector.infection_start_times[-1])
+            infectee.set_exposure_period(inf_to_exposed)
 
     def bind_population(self, population):
         super().bind_population(population)

pyEpiabm/pyEpiabm/tests/test_unit/test_core/test_person.py

@@ -149,6 +149,41 @@ def test_increment_secondary_infections(self):
         self.assertListEqual(self.person.secondary_infections_counts,
                              [2, 5, 2])
 
+    def test_set_exposure_period(self):
+        self.person.set_exposure_period(exposure_period=5)
+        self.assertEqual(self.person.exposure_period, 5)
+
+    def test_store_serial_interval_erroneous(self):
+        with self.assertRaises(RuntimeError) as ve_1:
+            self.person.store_serial_interval(4.0)
+        self.assertEqual(str(ve_1.exception),
+                         "Cannot call store_serial_interval while the"
+                         " exposure_period is None")
+        # Infect once with an exposure period for successful method call
+        self.person.set_exposure_period(1.0)
+        self.person.time_of_status_change = 2.0
+        self.person.store_serial_interval(latent_period=1.0)
+        # Do not add the exposure period for failure
+        self.person.time_of_status_change = 19.0
+        with self.assertRaises(RuntimeError) as ve_2:
+            self.person.store_serial_interval(latent_period=4.0)
+        self.assertEqual(str(ve_2.exception),
+                         "Cannot call store_serial_interval while the"
+                         " exposure_period is None")
+
+    def test_store_serial_interval(self):
+        self.person.set_exposure_period(1.0)
+        self.person.time_of_status_change = 2.0
+        self.person.store_serial_interval(latent_period=1.0)
+        self.person.set_exposure_period(2.0)
+        self.person.time_of_status_change = 6.0
+        self.person.store_serial_interval(latent_period=4.0)
+        self.person.set_exposure_period(5.0)
+        self.person.time_of_status_change = 19.0
+        self.person.store_serial_interval(latent_period=4.0)
+        self.assertDictEqual(self.person.serial_interval_dict,
+                             {0.0: [2.0, 6.0], 10.0: [9.0]})
+
 
 if __name__ == '__main__':
     unittest.main()