appointment_style_booking_models.qmd

---
title: Dealing With Appointment Bookings
code-annotations: hover
execute:
  eval: false
jupyter: python3
---

In many community-based services, it is important to model the process of booked appointments.

These are quite distinct from the processes we have modelled so far, where arrivals flow through the system immediately - or as quickly as possible, depending on the queues.

This works well for modelling a range of services, or smaller parts of more complex processes, such as

- emergency departments
- services without booking (e.g. a walk-in clinic)
- telephone helplines

However, in many services there is a process by which clients are booked into an appointment at some point in the future. There is often a delay of several days or weeks - which allows patients to receive communication about their appointment and make plans to attend it.

:::{.callout-note}
This video from NHS England explains why a certain level of waiting list can be a good thing for both patients and the service, and how to determine the ideal waiting list size for different services.

{{< video https://youtu.be/hga8HQFXaDw?si=KzNqZly5xf-bJ7S9>}}
:::

Appointment booking can take on additional layers of complexity, such as

- certain slots being set aside for more urgent referrals - and how to balance available capacity for these appointments with the importance of these clients being seen quickly
- an ongoing caseload of patients being held 'on the books' and returning for follow-up appointments at certain intervals
- a triage step prior to appointment booking where the referral may be deemed as inappropriate, giving another point at which referrals may leave the system (a **sink**)
- a certain proportion of patients not attending their appointment, with some exiting the system entirely while some may need to be rebooked

In this chapter, we will look at an implementation of a simple model of a mental health assessment service.

In this model, patients

- are referred to the service
- are booked in to the next available appointment
- wait until the appointment is carried out
- exit the system - we will make the assumption at this point that they are either discharged, or are referred on to a separate service for an intervention, but we will not model this part of the system

:::{.callout-note}
The contents of this section is based off the method employed by [Monks](https://orcid.org/0000-0003-2631-4481) in this [code](https://github.com/health-data-science-OR/stochastic_systems/tree/master/labs/simulation/lab5).
:::

## The appointment book

The key difference in this model is that we will feed in an additional object that represents the capacity of the clinic to see new clients.

To begin with, let's assume that

- there is a single clinic
- any client can be seen by any clinician
- they are open six days a week
- all appointments are the same length
- clients do not express any preference about being seen at a particular time of day
- everyone attends their appointment

```{python}
#| eval: true
#| echo: false

import pandas as pd
import simpy
import numpy as np
```

```{python}
#| eval: true
pd.read_csv("resources/shifts_simplest.csv")
```

Here, we have one row per day of the week. We will interpret an index of 0 as Monday and 6 as Sunday.

## Coding the example

### The g class

Rather than setting an interarrival time, we will instead set a value that represents the average annual demand for our clinic.

We will also pass in the dataframe of shifts.

Another new parameter is the minimum wait. To give patients time to receive their appointment letter and make a plan to attend, we don't want to just book the next available appointment, as this could be the very next day, with no time for clients to find out they are meant to be attending!

In the final new parameter, we will

Next, we return to parameters we have used before - the sim duration, which we have this time set as two years (365 days times 2). Note that compared to our previous model, where we interpreted each simpy time unit as 1 minute, we are now interpreting a single time unit as one *day*. We do not change anything in simpy itself to do this - but we just need to be careful that we remain consistent in our application of this throughout the rest of the model.

```{python}
#| eval: true

shifts_df = pd.read_csv("resources/shifts_simplest.csv")

# Class to store global parameter values.  We don't create an instance of this
# class - we just refer to the class blueprint itself to access the numbers
# inside.
class g:
    annual_demand = 3500
    shifts = shifts_df

    min_wait = 7

    sim_duration = 365 * 2
    number_of_runs = 100

```

### The Patient (entity) class

In our patient class, we will record an id as before.

We have a new attribute named 'booker' that we will create shortly.

We will also create a space to record the time patients arrive into the model, and the time they have their appointment.

```{python}
#| eval: true

# Class representing patients coming in to the clinic.
class Patient:
    def __init__(self, p_id, booker):
        self.id = p_id
        self.booker = booker

        self.arrival_time = 0
        self.waiting_time = 0
```


### The model class

#### The __init__method

We now need to make some important changes to the __init__method of the model, as well as create a few extra methods we can call on.

One of the key things we need to do is create two new dataframes based on our shift data (the dataframe of available daily slots). We have only provided the required information for a single week - our model will need to

- extrapolate this out into an array that covers the whole model (with a bit extra for appointments that are booked while the model is running, but are booked in for after the model has finished running)
- create a second array with the same dimensions that will be used to track the number of patients that have been booked in on a given day, allowing us to calculate if there are any slots still available when

:::{.callout-tip}
Note that we are using numpy throughout for the operations relating to the bookings. This is just due to the speed advantage of numpy in this context.
:::

When setting up our model class, we also want to create a distribution object that can be used to sample the number of daily arrivals to our system, based on the average number of yearly arrivals passed in our g class.

For more information on these functions, refer to #sec-reproducibility and #sec-distributions.

We will be using the `Poisson` class from the `sim-tools` package, so first need to run the following line.

```{python}
#| eval: true

from sim_tools.distributions import Poisson
```

```{python}
#| eval: true
class Model:
    # Constructor to set up the model for a run.  We pass in a run number when
    # we create a new model.
    def __init__(self, run_number):
        # Create a SimPy environment in which everything will live
        self.env = simpy.Environment()

        # Create a patient counter (which we'll use as a patient ID)
        self.patient_counter = 0

        # Store the passed in run number
        self.run_number = run_number

        ## NEW
        self.available_slots = None
        self.bookings = None

        ## NEW - run the new methods we have created below
        self.create_slots() ##NEW
        self.create_bookings() ##NEW

        ## NEW
        self.arrival_dist = Poisson(g.annual_demand / 52 / 7, # <1>
                                    random_seed=run_number*42) # <2>

        # Create a new Pandas DataFrame that will store some results against
        # the patient ID (which we'll use as the index).
        self.results_df = pd.DataFrame()
        self.results_df["Patient ID"] = [1]
        self.results_df["Q Time Appointment"] = [0.0] ##NEW
        self.results_df.set_index("Patient ID", inplace=True)

        # Create an attribute to store the mean waiting time for an appointment
        # across this run of the model
        self.mean_wait_time_appointment = 0 ## NEW
        self.mean_yearly_arrivals = 0 ## NEW

    ########################################
    ## ---------- NEW ------------------- ##
    ########################################

    def create_slots(self):

        available_slots = g.shifts.astype(np.uint8)
        template = available_slots.copy()

        #longer than run length as patients will need to book ahead
        for day in range(int((g.sim_duration/len(g.shifts))*3)):
            available_slots = pd.concat([available_slots, template.copy()],
                                         ignore_index=True)

        available_slots.index.rename('day', inplace=True)
        self.available_slots = available_slots

    def create_bookings(self):
        bookings = np.zeros(shape=(len(g.shifts), len(g.shifts.columns)), dtype=np.uint8)

        columns = [f'clinic_{i}' for i in range(1, len(g.shifts.columns)+1)]
        bookings_template = pd.DataFrame(bookings, columns=columns)

        bookings = bookings_template.copy()

        #longer than run length as patients will need to book ahead
        for day in range(int((g.sim_duration/len(g.shifts))*3)):
            bookings = pd.concat([bookings, bookings_template.copy()],
                                 ignore_index=True)

        bookings.index.rename('day', inplace=True)
        self.bookings = bookings

    ########################################
    ## ---------- END NEW --------------- ##
    ########################################
```
1. The parameter we pass to the poisson distribution is the average number of daily arrivals; this will be the average number of yearly arrivals divided by the number of weeks in a year and the number of days in a week.
2. This will ensure we have a different random pattern of arrivals per run, but that the number of arrivals is reproducible across trials.

Let's look at the outputs from this.

```{python}
#| eval: true
model = Model(run_number=1)

model.create_slots()
model.available_slots
```

```{python}
#| eval: true
model.create_bookings()
model.bookings
```

### The booker class

Before we continue making changes to our model class, we want to introduce a new class that will deal with patient bookings.

:::{.callout-tip}
While at this stage these methods could easily be incorporated elsewhere - such as into the model class itself - separating it out into its own class will give us more flexibility in the future when we wish to add in features such as pooling of clinic resources or different booking protocols for high and low priority patients.

In this case, this is also the motivation for adding in a 'clinic ID' parameter that is passed in when making the booking. While we only have a single clinic in this version of the model, allowing for a situation where we have multiple clinics a client could attend, and can make a choice to send them to whichever of these potential clinics have the earliest available appointment.
:::


```{python}
#| eval: true
class Booker():
    '''
    Booking class.
    '''
    def __init__(self, model):
        self.priority = 1
        self.model = model

    def find_slot(self, t, clinic_id):
        '''
        Finds a slot in a diary of available slot

        Params:
        ------
        t: int,
            current simulation time in days
            required to prevent booking an appointment
            in the past

        clinic_id: int,
            index of clinic to look up slots for

        Returns:
        -------
        (int, int)
        (booking_t, best_clinic_idx)

        '''
        # to reduce runtime - drop down to numpy
        available_slots_np = self.model.available_slots.to_numpy()

        # get the clinic slots t + min_wait forward
        clinic_slots = available_slots_np[t + g.min_wait: , clinic_id]

        # get the earliest day number (its the name of the series)
        best_t = np.where((clinic_slots.reshape(len(clinic_slots),1).sum(axis=1) > 0))[0][0]

        # Note that to get the index (day) of the actual booking, we
        # need to add the simulation time (t) and the minimum wait to the
        # index of the best time we found
        booking_t = t + g.min_wait + best_t

        return booking_t, clinic_id


    def book_slot(self, booking_t, clinic_id):
        '''
        Book a slot on day t for clinic c

        A slot is removed from args.available_slots
        A appointment is recorded in args.bookings.iat

        Params:
        ------
        booking_t: int
            Day of booking

        clinic_id: int
            the clinic identifier
        '''

        # Reduce the number of available slots by one at the point of the booking
        self.model.available_slots.iat[booking_t, clinic_id] -= 1

        # Increase the number of bookings we have on that day
        self.model.bookings.iat[booking_t, clinic_id] += 1

```


Let's take a look at the output of these.

```{python}
#| eval: true

sample_booker = Booker(model)

booking_t, clinic_id = sample_booker.find_slot(t=10, clinic_id=0)
print(f"Booking t: {booking_t}")
print(f"Clinic Index: {clinic_id}")
```

```{python}
#| eval: true

booking_t, clinic_id = sample_booker.find_slot(t=320, clinic_id=0)
print(f"Booking t: {booking_t}")
print(f"Clinic Index: {clinic_id}")
```


```{python}
#| eval: true

booking_t, clinic_id = sample_booker.find_slot(t=0, clinic_id=0)
print(f"Booking t: {booking_t}")
print(f"Clinic Index: {clinic_id}")

sample_booker.book_slot(booking_t=booking_t, clinic_id=clinic_id)

model.available_slots
```

```{python}
#| eval: true
model.bookings
```


### Further changes to the Model class

Now we have our booker method, we can make the remaining changes required to the model class.

#### The generator_patient_arrivals method

This method changes quite significantly from our existing models.

Instead of generating a patient, calculating the length of time that elapses until the next patient arrives and then waiting for that time to elapse, we mode through the simulation one day at a time, performing the following steps:

- calculating (sampling) the number of arrivals per day
- looping through each referral and creating a new patient object
- creating an instance of the booker class for each patient
- starting a referral process for that patient

When this is complete for every patient who is generated for that day, we can step forward to the next day by waiting for one time unit to elapse.

```{python}
def generator_patient_arrivals(self):
    for t in itertools.count():

        #total number of referrals today
        n_referrals = self.arrival_dist.sample()

        #loop through all referrals recieved that day
        for i in range(n_referrals):
            self.patient_counter += 1

            booker = Booker(model=self)

            # Create instance of Patient
            p = Patient(p_id=self.patient_counter, booker=booker)

            # Start a referral assessment process for patient.
            self.env.process(self.attend_clinic(p))

        #timestep by one day
        yield self.env.timeout(1)
```

#### The attend_clinic method

This method also needs to change quite significantly.

```{python}
def attend_clinic(self, patient):
    patient.arrival_time = self.env.now

    best_t, clinic_id = (
            patient.booker.find_slot(
              t = patient.arrival_time, # <1>
              clinic_id = 0 # <2>
              )
    )

    #book slot at clinic = time of referral + waiting_time
    #
    patient.booker.book_slot(best_t, clinic_id)

    # Wait for the time until the appointment to elapse
    yield self.env.timeout(best_t - patient.arrival_time)

    patient.waiting_time = self.env.now - patient.arrival_time

    self.results_df.at[patient.id, "Q Time Appointment"] = (
               patient.waiting_time
               )
```
1. The arrival time of the patient is passed in so that only appointments in the future are considered as eligible slots.
2. In this example we are only looking at a single clinic. In later versions of this model, patients who arrive will also have a preferred 'home' clinic, which will mean that different clinic IDs can be passed in at this stage.

#### The calculate_run_results method

Here, we are just altering the column name we refer to when calculating our metric of interest.

```{python}
    def calculate_run_results(self):
        # Take the mean of the queuing times for an appointment across this run of the model
        self.mean_wait_time_appointment= self.results_df["Q Time Appointment"].mean()
        # use our patient counter to track how many patients turn up on average during each
        # year of the simulation.
        self.mean_yearly_arrivals = self.patient_counter / (g.sim_duration / 365)
```

#### The run method

The run method is unchanged.

### The trial class

Our trial class is fundamentally unchanged - the main differences relate to the changes to the metrics we are interested in tracking.

```{python}
# Class representing a Trial for our simulation - a batch of simulation runs.
class Trial:
    # The constructor sets up a pandas dataframe that will store the key
    # results from each run (just the mean queuing time for the nurse here)
    # against run number, with run number as the index.
    def  __init__(self):
        self.df_trial_results = pd.DataFrame()
        self.df_trial_results["Run Number"] = [0]
        self.df_trial_results["Mean Appointment Wait (Days)"] = [0.0] ##NEW
        self.df_trial_results["Average Yearly Arrivals"] = [0.0] ##NEW
        self.df_trial_results.set_index("Run Number", inplace=True)

    # Method to print out the results from the trial.  In real world models,
    # you'd likely save them as well as (or instead of) printing them
    def print_trial_results(self):
        print ("Trial Results")
        print (self.df_trial_results)

    # Method to run a trial
    def run_trial(self):
        # Run the simulation for the number of runs specified in g class.
        # For each run, we create a new instance of the Model class and call its
        # run method, which sets everything else in motion.  Once the run has
        # completed, we grab out the stored run results (just mean queuing time
        # here) and store it against the run number in the trial results
        # dataframe.
        for run in range(g.number_of_runs):
            my_model = Model(run)
            my_model.run()

            self.df_trial_results.loc[run, "Mean Appointment Wait (Days)"] = [my_model.mean_wait_time_appointment] ##NEW
            self.df_trial_results.loc[run, "Average Yearly Arrivals"] = [my_model.mean_yearly_arrivals] ##NEW


        # Once the trial (ie all runs) has completed, print the final results
        self.print_trial_results()
```

### The full code

:::{.callout-note collapse="true"}
### Click here to view the full code


```{python}
#| eval: true

from sim_tools.distributions import Poisson
import pandas as pd
import numpy as np
import simpy
import itertools

shifts_df = pd.read_csv("resources/shifts_simplest.csv")

# Class to store global parameter values.  We don't create an instance of this
# class - we just refer to the class blueprint itself to access the numbers
# inside.
class g:
    annual_demand = 3500
    shifts = shifts_df

    min_wait = 7

    sim_duration = 365 * 2
    number_of_runs = 10

# Class representing patients coming in to the clinic.
class Patient:
    def __init__(self, p_id, booker):
        self.id = p_id
        self.booker = booker

        self.arrival_time = 0
        self.waiting_time = 0

class Booker():
    '''
    Booking class.
    '''
    def __init__(self, model):
        self.priority = 1
        self.model = model

    def find_slot(self, t, clinic_id):
        '''
        Finds a slot in a diary of available slot

        Params:
        ------
        t: int,
            current simulation time in days
            required to prevent booking an appointment
            in the past

        clinic_id: int,
            index of clinic to look up slots for

        Returns:
        -------
        (int, int)
        (booking_t, best_clinic_idx)

        '''
        # to reduce runtime - drop down to numpy
        available_slots_np = self.model.available_slots.to_numpy()

        # get the clinic slots t + min_wait forward
        clinic_slots = available_slots_np[t + g.min_wait: , clinic_id]

        # get the earliest day number (its the name of the series)
        best_t = np.where((clinic_slots.reshape(len(clinic_slots),1).sum(axis=1) > 0))[0][0]

        # Note that to get the index (day) of the actual booking, we
        # need to add the simulation time (t) and the minimum wait to the
        # index of the best time we found
        booking_t = t + g.min_wait + best_t

        return booking_t, clinic_id


    def book_slot(self, booking_t, clinic_id):
        '''
        Book a slot on day t for clinic c

        A slot is removed from args.available_slots
        A appointment is recorded in args.bookings.iat

        Params:
        ------
        booking_t: int
            Day of booking

        clinic_id: int
            the clinic identifier
        '''

        # Reduce the number of available slots by one at the point of the booking
        self.model.available_slots.iat[booking_t, clinic_id] -= 1

        # Increase the number of bookings we have on that day
        self.model.bookings.iat[booking_t, clinic_id] += 1

class Model:
    # Constructor to set up the model for a run.  We pass in a run number when
    # we create a new model.
    def __init__(self, run_number):
        # Create a SimPy environment in which everything will live
        self.env = simpy.Environment()

        # Create a patient counter (which we'll use as a patient ID)
        self.patient_counter = 0

        # Store the passed in run number
        self.run_number = run_number

        ## NEW
        self.available_slots = None
        self.bookings = None

        ## Populate these two items
        self.create_slots() ##NEW
        self.create_bookings() ##NEW


        ## NEW
        self.arrival_dist = Poisson(g.annual_demand / 52 / 7, # <1>
                                    random_seed=run_number*42) # <2>

        # Create a new Pandas DataFrame that will store some results against
        # the patient ID (which we'll use as the index).
        self.results_df = pd.DataFrame()
        self.results_df["Patient ID"] = [1]
        self.results_df["Q Time Appointment"] = [0.0] ##NEW
        self.results_df.set_index("Patient ID", inplace=True)

        # Create an attribute to store the mean waiting time for an appointment
        # across this run of the model
        self.mean_wait_time_appointment = 0 ## NEW
        self.mean_yearly_arrivals = 0 ## NEW

    ########################################
    ## ---------- NEW ------------------- ##
    ########################################

    def create_slots(self):

        available_slots = g.shifts.astype(np.uint8)
        template = available_slots.copy()

        #longer than run length as patients will need to book ahead
        for day in range(int((g.sim_duration/len(g.shifts))*3)):
            available_slots = pd.concat([available_slots, template.copy()],
                                         ignore_index=True)

        available_slots.index.rename('day', inplace=True)
        self.available_slots = available_slots

    def create_bookings(self):
        bookings = np.zeros(shape=(len(g.shifts), len(g.shifts.columns)), dtype=np.uint8)

        columns = [f'clinic_{i}' for i in range(1, len(g.shifts.columns)+1)]
        bookings_template = pd.DataFrame(bookings, columns=columns)

        bookings = bookings_template.copy()

        #longer than run length as patients will need to book ahead
        for day in range(int((g.sim_duration/len(g.shifts))*3)):
            bookings = pd.concat([bookings, bookings_template.copy()],
                                 ignore_index=True)

        bookings.index.rename('day', inplace=True)
        self.bookings = bookings

    def generator_patient_arrivals(self):
        for t in itertools.count():

            #total number of referrals today
            n_referrals = self.arrival_dist.sample()

            #loop through all referrals recieved that day
            for i in range(n_referrals):
                self.patient_counter += 1

                booker = Booker(model=self)

                # Create instance of Patient
                p = Patient(p_id=self.patient_counter, booker=booker)

                # Start a referral assessment process for patient.
                self.env.process(self.attend_clinic(p))

            #timestep by one day
            yield self.env.timeout(1)

    def attend_clinic(self, patient):
        patient.arrival_time = self.env.now

        best_t, clinic_id = (
                patient.booker.find_slot(
                  t = patient.arrival_time, # <1>
                  clinic_id = 0 # <2>
                  )
        )

        #book slot at clinic = time of referral + waiting_time
        patient.booker.book_slot(best_t, clinic_id)

        # Wait for the time until the appointment to elapse
        yield self.env.timeout(best_t - patient.arrival_time)

        patient.waiting_time = self.env.now - patient.arrival_time

        self.results_df.at[patient.id, "Q Time Appointment"] = (
                  patient.waiting_time
                  )

    def calculate_run_results(self):
        # Take the mean of the queuing times for the nurse across patients in
        # this run of the model.
        self.mean_wait_time_appointment= self.results_df["Q Time Appointment"].mean()
        self.mean_yearly_arrivals = self.patient_counter / (g.sim_duration / 365)

    # The run method starts up the DES entity generators, runs the simulation,
    # and in turns calls anything we need to generate results for the run
    def run(self):

        # Start up our DES entity generators that create new patients.  We've
        # only got one in this model, but we'd need to do this for each one if
        # we had multiple generators.
        self.env.process(self.generator_patient_arrivals())

        # Run the model for the duration specified in g class
        self.env.run(until=g.sim_duration)

        # Now the simulation run has finished, call the method that calculates
        # run results
        self.calculate_run_results()

        # Print the run number with the patient-level results from this run of
        # the model
        # Commented out for now
        #print (f"Run Number {self.run_number}")
        #print (self.results_df)

    ########################################
    ## ---------- END NEW --------------- ##
    ########################################

# Class representing a Trial for our simulation - a batch of simulation runs.
class Trial:
    # The constructor sets up a pandas dataframe that will store the key
    # results from each run (just the mean queuing time for the nurse here)
    # against run number, with run number as the index.
    def  __init__(self):
        self.df_trial_results = pd.DataFrame()
        self.df_trial_results["Run Number"] = [0]
        self.df_trial_results["Mean Appointment Wait (Days)"] = [0.0] ##NEW
        self.df_trial_results["Average Yearly Arrivals"] = [0.0] ##NEW
        self.df_trial_results.set_index("Run Number", inplace=True)

    # Method to print out the results from the trial.  In real world models,
    # you'd likely save them as well as (or instead of) printing them
    def print_trial_results(self):
        print ("Trial Results")
        print (self.df_trial_results)

    # Method to run a trial
    def run_trial(self):
        # Run the simulation for the number of runs specified in g class.
        # For each run, we create a new instance of the Model class and call its
        # run method, which sets everything else in motion.  Once the run has
        # completed, we grab out the stored run results (just mean queuing time
        # here) and store it against the run number in the trial results
        # dataframe.
        for run in range(g.number_of_runs):
            my_model = Model(run)
            my_model.run()

            self.df_trial_results.loc[run, "Mean Appointment Wait (Days)"] = [my_model.mean_wait_time_appointment] ##NEW
            self.df_trial_results.loc[run, "Average Yearly Arrivals"] = [my_model.mean_yearly_arrivals] ##NEW

        # Once the trial (ie all runs) has completed, print the final results
        self.print_trial_results()
```
:::

### Evaluating the results

```{python}
#| eval: true

# Create an instance of the Trial class
my_trial = Trial()

# Call the run_trial method of our Trial object
my_trial.run_trial()
```

## Tracking additional metrics

We may find it useful to understand how many of our available appointment slots are going unused in the simulation.

To do this, we can slice our available_slots and bookings objects to just include the period of interest.

Remember that available_slots refers to the number of remaining slots after bookings have been made - not the total available theoretical slots per day, which is stored in our g class as `shifts` - but remember that the shifts object is only a template for a seven day period rather than encompassing the whole model duration.

Therefore, to get the total number of possible slots, we must do `available_slots + bookings`.

This will add up the relevant values on a day-by-day basis.

We can then sum `available_slots` to get the total number of slots that weren't utilised, and then sum the result of `available_slots + bookings` to get the total number of slots that were available. By dividing the first result by the second, we get an indication of what proportion of slots were actually used.

:::{.callout-note}
Changes to the code are marked with ##NEW below
:::

```{python}
#| eval: true

class Model:
    # Constructor to set up the model for a run.  We pass in a run number when
    # we create a new model.
    def __init__(self, run_number):
        # Create a SimPy environment in which everything will live
        self.env = simpy.Environment()

        # Create a patient counter (which we'll use as a patient ID)
        self.patient_counter = 0

        # Store the passed in run number
        self.run_number = run_number

        self.available_slots = None
        self.bookings = None

        ## Populate these two items
        self.create_slots()
        self.create_bookings()

        self.arrival_dist = Poisson(g.annual_demand / 52 / 7, # <1>
                                    random_seed=run_number*42) # <2>

        # Create a new Pandas DataFrame that will store some results against
        # the patient ID (which we'll use as the index).
        self.results_df = pd.DataFrame()
        self.results_df["Patient ID"] = [1]
        self.results_df["Q Time Appointment"] = [0.0]
        self.results_df.set_index("Patient ID", inplace=True)

        # Create an attribute to store the mean waiting time for an appointment
        # across this run of the model
        self.mean_wait_time_appointment = 0
        self.mean_yearly_arrivals = 0
        self.percentage_slots_used = 0.0 ##NEW

    def create_slots(self):

        available_slots = g.shifts.astype(np.uint8)
        template = available_slots.copy()

        #longer than run length as patients will need to book ahead
        for day in range(int((g.sim_duration/len(g.shifts))*3)):
            available_slots = pd.concat([available_slots, template.copy()],
                                         ignore_index=True)

        available_slots.index.rename('day', inplace=True)
        self.available_slots = available_slots

    def create_bookings(self):
        bookings = np.zeros(shape=(len(g.shifts), len(g.shifts.columns)), dtype=np.uint8)

        columns = [f'clinic_{i}' for i in range(1, len(g.shifts.columns)+1)]
        bookings_template = pd.DataFrame(bookings, columns=columns)

        bookings = bookings_template.copy()

        #longer than run length as patients will need to book ahead
        for day in range(int((g.sim_duration/len(g.shifts))*3)):
            bookings = pd.concat([bookings, bookings_template.copy()],
                                 ignore_index=True)

        bookings.index.rename('day', inplace=True)
        self.bookings = bookings

    def generator_patient_arrivals(self):
        for t in itertools.count():

            #total number of referrals today
            n_referrals = self.arrival_dist.sample()

            #loop through all referrals recieved that day
            for i in range(n_referrals):
                self.patient_counter += 1

                booker = Booker(model=self)

                # Create instance of Patient
                p = Patient(p_id=self.patient_counter, booker=booker)

                # Start a referral assessment process for patient.
                self.env.process(self.attend_clinic(p))

            #timestep by one day
            yield self.env.timeout(1)

    def attend_clinic(self, patient):
        patient.arrival_time = self.env.now

        best_t, clinic_id = (
                patient.booker.find_slot(
                  t = patient.arrival_time, # <1>
                  clinic_id = 0 # <2>
                  )
        )

        #book slot at clinic = time of referral + waiting_time
        patient.booker.book_slot(best_t, clinic_id)

        # Wait for the time until the appointment to elapse
        yield self.env.timeout(best_t - patient.arrival_time)

        patient.waiting_time = self.env.now - patient.arrival_time

        self.results_df.at[patient.id, "Q Time Appointment"] = (
                  patient.waiting_time
                  )

    def calculate_run_results(self):
        # Take the mean of the queuing times for the nurse across patients in
        # this run of the model.
        self.mean_wait_time_appointment= self.results_df["Q Time Appointment"].mean()
        self.mean_yearly_arrivals = self.patient_counter / (g.sim_duration / 365)
        #########
        ## NEW ##
        #########

        slots_unused = self.available_slots.clinic_1.values[ : g.sim_duration]

        slots_used = self.bookings.clinic_1.values[ : g.sim_duration]

        total_slots_available_daily = np.add(slots_unused, slots_used)

        self.percentage_slots_used = (sum(slots_used) / sum(total_slots_available_daily))

    # The run method starts up the DES entity generators, runs the simulation,
    # and in turns calls anything we need to generate results for the run
    def run(self):

        # Start up our DES entity generators that create new patients.  We've
        # only got one in this model, but we'd need to do this for each one if
        # we had multiple generators.
        self.env.process(self.generator_patient_arrivals())

        # Run the model for the duration specified in g class
        self.env.run(until=g.sim_duration)

        # Now the simulation run has finished, call the method that calculates
        # run results
        self.calculate_run_results()

        # Print the run number with the patient-level results from this run of
        # the model
        # Commented out for now
        #print (f"Run Number {self.run_number}")
        #print (self.results_df)

class Trial:
    # The constructor sets up a pandas dataframe that will store the key
    # results from each run (just the mean queuing time for the nurse here)
    # against run number, with run number as the index.
    def  __init__(self):
        self.df_trial_results = pd.DataFrame()
        self.df_trial_results["Run Number"] = [0]
        self.df_trial_results["Mean Appointment Wait (Days)"] = [0.0] ##NEW
        self.df_trial_results["Average Yearly Arrivals"] = [0.0] ##NEW
        self.df_trial_results["Percentage of Slots Used"] = [0.0]
        self.df_trial_results.set_index("Run Number", inplace=True)

    # Method to print out the results from the trial.  In real world models,
    # you'd likely save them as well as (or instead of) printing them
    def print_trial_results(self):
        print ("Trial Results")
        print (self.df_trial_results)

    # Method to run a trial
    def run_trial(self):
        # Run the simulation for the number of runs specified in g class.
        # For each run, we create a new instance of the Model class and call its
        # run method, which sets everything else in motion.  Once the run has
        # completed, we grab out the stored run results (just mean queuing time
        # here) and store it against the run number in the trial results
        # dataframe.
        for run in range(g.number_of_runs):
            my_model = Model(run)
            my_model.run()

            self.df_trial_results.loc[run, "Mean Appointment Wait (Days)"] = [my_model.mean_wait_time_appointment]

            self.df_trial_results.loc[run, "Average Yearly Arrivals"] = [my_model.mean_yearly_arrivals]

            self.df_trial_results.loc[run, "Percentage of Slots Used"] = [my_model.percentage_slots_used] ##NEW

        # Once the trial (ie all runs) has completed, print the final results
        self.print_trial_results()
```

### Evaluating the Results

```{python}
#| eval: true

# Create an instance of the Trial class
my_trial = Trial()

# Call the run_trial method of our Trial object
my_trial.run_trial()
```

:::{.callout-important}
Note that in this example we have not made use of a warm-up period.

In the early phases of the model, we would expect the utilisation of available slots to be lower, as there will be no booked appointments at the beginning of the model, and the first people to be booked in will be booked in a week ahead due to the mandatory delay period. Therefore the utilisation is likely to be an underestimate.

To learn more about setting warm-up periods, refer to #sec-warmup
:::

## Evaluating the impact of changing the number of yearly attendances
To make it possible to see the impact of changing the number of attendances, we will modify our run_trial method to return the final dataframe instead of printing it.

```{python}
#| eval: true
#|
class Trial:
    # The constructor sets up a pandas dataframe that will store the key
    # results from each run (just the mean queuing time for the nurse here)
    # against run number, with run number as the index.
    def  __init__(self):
        self.df_trial_results = pd.DataFrame()
        self.df_trial_results["Run Number"] = [0]
        self.df_trial_results["Mean Appointment Wait (Days)"] = [0.0] ##NEW
        self.df_trial_results["Average Yearly Arrivals"] = [0.0] ##NEW
        self.df_trial_results["Percentage of Slots Used"] = [0.0]
        self.df_trial_results.set_index("Run Number", inplace=True)

    # Method to run a trial
    def run_trial(self):
        # Run the simulation for the number of runs specified in g class.
        # For each run, we create a new instance of the Model class and call its
        # run method, which sets everything else in motion.  Once the run has
        # completed, we grab out the stored run results (just mean queuing time
        # here) and store it against the run number in the trial results
        # dataframe.
        for run in range(g.number_of_runs):
            my_model = Model(run)
            my_model.run()

            self.df_trial_results.loc[run, "Mean Appointment Wait (Days)"] = [my_model.mean_wait_time_appointment]

            self.df_trial_results.loc[run, "Average Yearly Arrivals"] = [my_model.mean_yearly_arrivals]

            self.df_trial_results.loc[run, "Percentage of Slots Used"] = [my_model.percentage_slots_used] ##NEW

        # Once the trial (ie all runs) has completed, print the final results
        return self.df_trial_results
```

We can then iterate through and store the results for a number of different average yearly demand.

```{python}
#| eval: true
#|
yearly_demand = [2500, 2750, 3000, 3250, 3500, 3750, 4000]
results_frames = []

# We will also extend the duration of each run to 4 years
sim_duration = 365 * 4

for demand in yearly_demand:
    g.annual_demand = demand

    # Create an instance of the Trial class
    my_trial = Trial()

    # Call the run_trial method of our Trial object
    results_df = my_trial.run_trial()

    results_df['Demand'] = demand

    results_frames.append(results_df)

final_results_df = pd.concat(results_frames)
```

We can then visualise the results in bar charts.

First, to make it easier to visualise multiple metrics at once, we will reshape our results dataframe from wide to long.

The original dataframe looks like this:

```{python}
#| eval: true

final_results_df
```

After using the melt function of pandas, our dataframe now has one row per run.

```{python}
#| eval: true

final_results_df_long = pd.melt(
  final_results_df.reset_index(),
  id_vars=["Demand", "Run Number"]
  )

final_results_df_long
```

This allows us to use the `facet_rows` argument in plotly express.

```{python}
#| eval: true

import plotly.express as px

fig = px.box(
  final_results_df_long,
  x="Demand",
  y="value",
  facet_col="variable",
  facet_col_wrap=1,
  height=1200
)

fig.update_yaxes(matches=None)
fig.for_each_annotation(lambda a: a.update(text=a.text.split("=")[-1]))
fig.show()
```