!pip install simpy
!pip install sim_tools
import simpy
from sim_tools.distributions import Exponential, Lognormal
import pandas as pd
import math
from scipy import stats
import numpy as npPart A
In your groups, study each code cell. In the cell above it, double click and write a description of what you think that block of code is doing. Try to work out what each line of code is doing.
Make sure you run the following cell first - this will install and import the necessary packages into this session. You do not need to write an explanation of this code block.
The following code block is…
class Patient:
def __init__(self, p_id):
self.id = p_id
self.q_time_registration = pd.NA
self.q_time_nurse = pd.NAThe following code block is…
class Param:
def __init__(
self,
mean_patient_inter = 3,
mean_registration_time = 3,
sd_registration_time = 0.5,
mean_nurse_consult_time = 6,
sd_nurse_consult_time = 1,
num_receptionists = 1,
num_nurses = 1,
sim_duration = 120,
num_replications = 5
):
self.mean_patient_inter = mean_patient_inter
self.mean_registration_time = mean_registration_time
self.sd_registration_time = sd_registration_time
self.mean_nurse_consult_time = mean_nurse_consult_time
self.sd_nurse_consult_time = sd_nurse_consult_time
self.num_receptionists = num_receptionists
self.num_nurses = num_nurses
self.sim_duration = sim_duration
self.num_replications = num_replicationsIMPORTANT NOTE : The following code cells split the large Model class definition over multiple cells. Because of this, we use something in Object Oriented Programming known as Inheritance, which is where we say that a class is like another class but with these differences. It’s used here to add a new method to the same Model class in each code block - hence you will see
class Model(Model)as the first line on each code cell after the first one. You wouldn’t normally do this in a real world model - we’re only doing it here to ensure that the code cells all work independently. So you can ignore that first line in each and just interpret the rest of the code in that cell as another method in the Model class definition.
Later, you’ll see we do the same for the Trial class.
The following code block is…
class Model:
def __init__(self, param, replication_id):
self.param = param
self.replication_id = replication_id
self.env = simpy.Environment()
self.patient_counter = 0
self.receptionist = simpy.Resource(
self.env, capacity=self.param.num_receptionists
)
self.nurse = simpy.Resource(self.env, capacity=self.param.num_nurses)
ss = np.random.SeedSequence(self.replication_id)
seeds = ss.spawn(3)
self.patient_inter_dist = Exponential(
mean=self.param.mean_patient_inter,
random_seed=seeds[0]
)
self.registration_time_dist = Lognormal(
mean=self.param.mean_registration_time,
stdev=self.param.sd_registration_time,
random_seed=seeds[1]
)
self.nurse_consult_time_dist = Lognormal(
mean=self.param.mean_nurse_consult_time,
stdev=self.param.sd_nurse_consult_time,
random_seed=seeds[2]
)
self.list_of_patients = []
self.mean_q_time_registration = pd.NA
self.sd_q_time_registration = pd.NA
self.perc_90_q_time_registration = pd.NA
self.mean_q_time_nurse = pd.NA
self.sd_q_time_nurse = pd.NA
self.perc_90_q_time_nurse = pd.NAThe following code block is…
class Model(Model):
def attend_clinic(self, patient):
start_q_registration = self.env.now
with self.receptionist.request() as req:
yield req
end_q_registration = self.env.now
patient.q_time_registration = (
end_q_registration - start_q_registration
)
sampled_reg_act_time = self.registration_time_dist.sample()
yield self.env.timeout(sampled_reg_act_time)
start_q_nurse = self.env.now
with self.nurse.request() as req:
yield req
end_q_nurse = self.env.now
patient.q_time_nurse = end_q_nurse - start_q_nurse
sampled_nurse_act_time = self.nurse_consult_time_dist.sample()
yield self.env.timeout(sampled_nurse_act_time)The following code block is…
class Model(Model):
def generator_patient_arrivals(self):
while True:
self.patient_counter += 1
p = Patient(self.patient_counter)
self.list_of_patients.append(p)
self.env.process(self.attend_clinic(p))
sampled_inter = self.patient_inter_dist.sample()
yield self.env.timeout(sampled_inter)The following code block is…
class Model(Model):
def run_model(self):
self.env.process(self.generator_patient_arrivals())
self.env.run(until=self.param.sim_duration)I’ll give you a couple tips on this next one as it’s tough if you don’t know Python.
Tip 1 : In Python there is something know as list comprehension. This is where you can use a shortcut to create a list based on another list, potentially where certain criteria are met too. For example :
list_b = [patient.age for patient in list_a if age>60]would create a second list (list_b) that contains only patient ages (and only those over 60) from list_a
Tip 2 : object.__dict__ extracts all of the attributes belonging to an object
The following code block is…
class Model(Model):
def convert_entity_list_to_dataframe(self, entity_list):
entity_dateframe = pd.DataFrame(
entity.__dict__ for entity in entity_list
)
return entity_dateframeThe following code block is…
class Model(Model):
def calculate_run_results(self, entity_dataframe):
self.mean_q_time_registration = (
entity_dataframe["q_time_registration"].mean()
)
self.sd_q_time_registration = (
entity_dataframe["q_time_registration"].std()
)
self.perc_90_q_time_registration = (
entity_dataframe["q_time_registration"].quantile(0.9)
)
self.mean_q_time_nurse = (
entity_dataframe["q_time_nurse"].mean()
)
self.sd_q_time_nurse = entity_dataframe["q_time_nurse"].std()
self.perc_90_q_time_nurse = (
entity_dataframe["q_time_nurse"].quantile(0.9)
)The following code block is…
class Trial:
def __init__(self, param):
self.param = param
self.list_of_simulation_replications = []
self.trial_mean_q_time_registration = pd.NA
self.trial_sd_q_time_registration = pd.NA
self.trial_perc_90_q_time_registration = pd.NA
self.trial_mean_q_time_nurse = pd.NA
self.trial_sd_q_time_nurse = pd.NA
self.trial_perc_90_q_time_nurse = pd.NA
self.ci_lower_q_time_registration = pd.NA
self.ci_upper_q_time_registration = pd.NA
self.se_q_time_registration = pd.NA
self.ci_lower_q_time_nurse = pd.NA
self.ci_upper_q_time_nurse = pd.NA
self.se_q_time_nurse = pd.NAThe following code block is…
class Trial(Trial):
def run_trial(self):
for replication_id in range(self.param.num_replications):
model_replication = Model(self.param, replication_id)
model_replication.run_model()
patient_df = model_replication.convert_entity_list_to_dataframe(
model_replication.list_of_patients
)
model_replication.calculate_run_results(patient_df)
self.list_of_simulation_replications.append(model_replication)The following code block is…
class Trial(Trial):
def calculate_trial_results(self):
self.replication_df = pd.DataFrame(
replication.__dict__ for replication in
self.list_of_simulation_replications
)
self.trial_mean_q_time_registration = (
self.replication_df["mean_q_time_registration"].mean()
)
self.trial_sd_q_time_registration = (
self.replication_df["mean_q_time_registration"].std()
)
self.trial_perc_90_q_time_registration = (
self.replication_df["mean_q_time_registration"].quantile(0.9)
)
self.trial_mean_q_time_nurse = (
self.replication_df["mean_q_time_nurse"].mean()
)
self.trial_sd_q_time_nurse = (
self.replication_df["mean_q_time_nurse"].std()
)
self.trial_perc_90_q_time_nurse = (
self.replication_df["mean_q_time_nurse"].quantile(0.9)
)
self.se_q_time_registration = (
self.trial_sd_q_time_registration / math.sqrt(
self.param.num_replications
)
)
self.se_q_time_nurse = (
self.trial_sd_q_time_nurse / math.sqrt(self.param.num_replications)
)
t = stats.t.ppf(0.975, df=self.param.num_replications-1)
self.ci_lower_q_time_registration = (
self.trial_mean_q_time_registration - (
t * self.se_q_time_registration
)
)
self.ci_upper_q_time_registration = (
self.trial_mean_q_time_registration + (
t * self.se_q_time_registration
)
)
self.ci_lower_q_time_nurse = (
self.trial_mean_q_time_nurse - (t * self.se_q_time_nurse)
)
self.ci_upper_q_time_nurse = (
self.trial_mean_q_time_nurse + (t * self.se_q_time_nurse)
)The following code block is…
base_case_params = Param()
base_case_trial = Trial(base_case_params)
base_case_trial.run_trial()
base_case_trial.calculate_trial_results()
what_if_params = Param(num_nurses=2, num_receptionists=2)
what_if_trial = Trial(what_if_params)
what_if_trial.run_trial()
what_if_trial.calculate_trial_results()The following code block is…
print ("BASE CASE TRIAL RESULTS")
print ("-----------------------")
print ("Queuing Time for Registration")
print (f"Mean : {base_case_trial.trial_mean_q_time_registration:.2f} minutes")
print (f"SD : {base_case_trial.trial_sd_q_time_registration:.2f} minutes")
print (
f"90th Perc : {base_case_trial.trial_perc_90_q_time_registration:.2f}",
"minutes"
)
print (f"Standard Error : {base_case_trial.se_q_time_registration:.2f}")
print (
f"95% CI : ({base_case_trial.ci_lower_q_time_registration:.2f}, ",
f"{base_case_trial.ci_upper_q_time_registration:.2f}) minutes"
)
print ()
print ("Queuing Time for the Nurse")
print (f"Mean : {base_case_trial.trial_mean_q_time_nurse:.2f} minutes")
print (f"SD : {base_case_trial.trial_sd_q_time_nurse:.2f} minutes")
print (f"90th Perc : {base_case_trial.trial_perc_90_q_time_nurse:.2f} minutes")
print (f"Standard Error : {base_case_trial.se_q_time_nurse:.2f}")
print (
f"95% CI : ({base_case_trial.ci_lower_q_time_nurse:.2f}, ",
f"{base_case_trial.ci_upper_q_time_nurse:.2f}) minutes"
)
print ()
print ("2 NURSES, 2 RECEPTIONISTS TRIAL RESULTS")
print ("----------------------")
print ("Queuing Time for Registration")
print (f"Mean : {what_if_trial.trial_mean_q_time_registration:.2f} minutes")
print (f"SD : {what_if_trial.trial_sd_q_time_registration:.2f} minutes")
print (
f"90th Perc : {what_if_trial.trial_perc_90_q_time_registration:.2f}",
"minutes"
)
print (f"Standard Error : {what_if_trial.se_q_time_registration:.2f}")
print (
f"95% CI : ({what_if_trial.ci_lower_q_time_registration:.2f}, ",
f"{what_if_trial.ci_upper_q_time_registration:.2f}) minutes"
)
print ()
print ("Queuing Time for the Nurse")
print (f"Mean : {what_if_trial.trial_mean_q_time_nurse:.2f} minutes")
print (f"SD : {what_if_trial.trial_sd_q_time_nurse:.2f} minutes")
print (f"90th Perc : {what_if_trial.trial_perc_90_q_time_nurse:.2f} minutes")
print (f"Standard Error : {what_if_trial.se_q_time_nurse:.2f}")
print (
f"95% CI : ({what_if_trial.ci_lower_q_time_nurse:.2f}, ",
f"{what_if_trial.ci_upper_q_time_nurse:.2f}) minutes"
)
print ()Part B
After you have interpreted and run all of the code cells above, try the following tasks.
- Write a second “What If?” scenario that tests having 1 receptionist, 2 nurses, but average (and standard deviation) of registration and nurse consultation times tripled. Run this scenario, then print the results.
# INSERT CODE FOR TASK 1 BELOW- Change the pathway being modelled so that, after seeing the nurse, there is a 40% chance that a patient will go on to see a doctor. Patients who don’t simply leave. The mean consultation time with the doctor is 10 minutes, and has a standard deviation of 3 minutes. There is currently 1 doctor.
In each of the code cells below, insert or modify code where indicated, and make sure you run the cells after making changes.
Here are some things that will be useful for you to know for this activity :
- It is good practice not to hard code things like probabilities in your models. You should consider making it a changeable parameter.
- We can use conditional logic to determine whether a block of indented code runs or not. Here’s an example of a piece of Python code that would print a different message depending on if the patient’s age is less than 18 or not :
if patient.age < 18:
print ("You're young!")
else:
print ("You're young at heart, I'm sure")- If statements do not have to include an else condition
- One way of deciding something randomly is to sample a random number between 0 and 1. If the number sampled is less than the probability of the event happening, then we say it happened. This works because the probability of picking a value between 0 and the given probability value is exactly equal to the probability value when sampling from a uniform distribution (where all possibilities are equally likely)
- If we want our models to be reproducible, then EVERY point of randomness in our models must be seeded with a random seed
- The NumPy package contains something known as a default_rng, which is a random number generator (RNG) from which we can sample numbers, and provide a random seed. The following two lines of code will create a new default_rng (line 1), and sample a value between 0 and 1 from that RNG (line 2) :
name_of_my_rng = np.random.default_rng(my_random_seed)
sampled_value = name_of_my_rng.random()class Patient:
def __init__(self, p_id):
self.id = p_id
self.q_time_registration = pd.NA
self.q_time_nurse = pd.NA
# WRITE YOUR NEW CODE BELOWclass Param:
def __init__(
self,
mean_patient_inter = 5,
mean_registration_time = 3,
sd_registration_time = 0.5,
mean_nurse_consult_time = 6,
sd_nurse_consult_time = 1,
num_receptionists = 1,
num_nurses = 1,
sim_duration = 120,
num_replications = 5,
# WRITE YOUR NEW CODE BELOW
):
self.mean_patient_inter = mean_patient_inter
self.mean_registration_time = mean_registration_time
self.sd_registration_time = sd_registration_time
self.mean_nurse_consult_time = mean_nurse_consult_time
self.sd_nurse_consult_time = sd_nurse_consult_time
self.num_receptionists = num_receptionists
self.num_nurses = num_nurses
self.sim_duration = sim_duration
self.num_replications = num_replications
# WRITE YOUR NEW CODE BELOWclass Model:
def __init__(self, param, replication_id):
self.param = param
self.replication_id = replication_id
self.env = simpy.Environment()
self.patient_counter = 0
self.receptionist = simpy.Resource(
self.env, capacity=self.param.num_receptionists
)
self.nurse = simpy.Resource(self.env, capacity=self.param.num_nurses)
# WRITE YOUR NEW CODE BELOW
ss = np.random.SeedSequence(self.replication_id)
# ADAPT THE LINE OF CODE BELOW
seeds = ss.spawn(3)
self.patient_inter_dist = Exponential(
mean=self.param.mean_patient_inter,
random_seed=seeds[0]
)
self.registration_time_dist = Lognormal(
mean=self.param.mean_registration_time,
stdev=self.param.sd_registration_time,
random_seed=seeds[1]
)
self.nurse_consult_time_dist = Lognormal(
mean=self.param.mean_nurse_consult_time,
stdev=self.param.sd_nurse_consult_time,
random_seed=seeds[2]
)
# WRITE YOUR NEW CODE BELOW
self.list_of_patients = []
self.mean_q_time_registration = pd.NA
self.sd_q_time_registration = pd.NA
self.perc_90_q_time_registration = pd.NA
self.mean_q_time_nurse = pd.NA
self.sd_q_time_nurse = pd.NA
self.perc_90_q_time_nurse = pd.NA
# WRITE YOUR NEW CODE BELOW
def attend_clinic(self, patient):
start_q_registration = self.env.now
with self.receptionist.request() as req:
yield req
end_q_registration = self.env.now
patient.q_time_registration = (
end_q_registration - start_q_registration
)
sampled_reg_act_time = self.registration_time_dist.sample()
yield self.env.timeout(sampled_reg_act_time)
start_q_nurse = self.env.now
with self.nurse.request() as req:
yield req
end_q_nurse = self.env.now
patient.q_time_nurse = end_q_nurse - start_q_nurse
sampled_nurse_act_time = self.nurse_consult_time_dist.sample()
yield self.env.timeout(sampled_nurse_act_time)
# WRITE YOUR NEW CODE BELOW
def generator_patient_arrivals(self):
while True:
self.patient_counter += 1
p = Patient(self.patient_counter)
self.list_of_patients.append(p)
self.env.process(self.attend_clinic(p))
sampled_inter = self.patient_inter_dist.sample()
yield self.env.timeout(sampled_inter)
def run_model(self):
self.env.process(self.generator_patient_arrivals())
self.env.run(until=self.param.sim_duration)
def convert_entity_list_to_dataframe(self, entity_list):
entity_dateframe = pd.DataFrame(
entity.__dict__ for entity in entity_list
)
return entity_dateframe
def calculate_run_results(self, entity_dataframe):
self.mean_q_time_registration = (
entity_dataframe["q_time_registration"].mean()
)
self.sd_q_time_registration = (
entity_dataframe["q_time_registration"].std()
)
self.perc_90_q_time_registration = (
entity_dataframe["q_time_registration"].quantile(0.9)
)
self.mean_q_time_nurse = (
entity_dataframe["q_time_nurse"].mean()
)
self.sd_q_time_nurse = entity_dataframe["q_time_nurse"].std()
self.perc_90_q_time_nurse = (
entity_dataframe["q_time_nurse"].quantile(0.9)
)
# WRITE YOUR NEW CODE BELOWclass Trial:
def __init__(self, param):
self.param = param
self.list_of_simulation_replications = []
self.trial_mean_q_time_registration = pd.NA
self.trial_sd_q_time_registration = pd.NA
self.trial_perc_90_q_time_registration = pd.NA
self.trial_mean_q_time_nurse = pd.NA
self.trial_sd_q_time_nurse = pd.NA
self.trial_perc_90_q_time_nurse = pd.NA
self.ci_lower_q_time_registration = pd.NA
self.ci_upper_q_time_registration = pd.NA
self.se_q_time_registration = pd.NA
self.ci_lower_q_time_nurse = pd.NA
self.ci_upper_q_time_nurse = pd.NA
self.se_q_time_nurse = pd.NA
# WRITE YOUR NEW CODE BELOW
def run_trial(self):
for replication_id in range(self.param.num_replications):
model_replication = Model(self.param, replication_id)
model_replication.run_model()
patient_df = model_replication.convert_entity_list_to_dataframe(
model_replication.list_of_patients
)
model_replication.calculate_run_results(patient_df)
self.list_of_simulation_replications.append(model_replication)
def calculate_trial_results(self):
self.replication_df = pd.DataFrame(
replication.__dict__ for replication in
self.list_of_simulation_replications
)
self.trial_mean_q_time_registration = (
self.replication_df["mean_q_time_registration"].mean()
)
self.trial_sd_q_time_registration = (
self.replication_df["mean_q_time_registration"].std()
)
self.trial_perc_90_q_time_registration = (
self.replication_df["mean_q_time_registration"].quantile(0.9)
)
self.trial_mean_q_time_nurse = (
self.replication_df["mean_q_time_nurse"].mean()
)
self.trial_sd_q_time_nurse = (
self.replication_df["mean_q_time_nurse"].std()
)
self.trial_perc_90_q_time_nurse = (
self.replication_df["mean_q_time_nurse"].quantile(0.9)
)
# WRITE YOUR NEW CODE BELOW
self.se_q_time_registration = (
self.trial_sd_q_time_registration / math.sqrt(
self.param.num_replications
)
)
self.se_q_time_nurse = (
self.trial_sd_q_time_nurse / math.sqrt(self.param.num_replications)
)
# WRITE YOUR NEW CODE BELOW
t = stats.t.ppf(0.975, df=self.param.num_replications-1)
self.ci_lower_q_time_registration = (
self.trial_mean_q_time_registration - (
t * self.se_q_time_registration
)
)
self.ci_upper_q_time_registration = (
self.trial_mean_q_time_registration + (
t * self.se_q_time_registration
)
)
self.ci_lower_q_time_nurse = (
self.trial_mean_q_time_nurse - (t * self.se_q_time_nurse)
)
self.ci_upper_q_time_nurse = (
self.trial_mean_q_time_nurse + (t * self.se_q_time_nurse)
)
# WRITE YOUR NEW CODE BELOWbase_case_params = Param()
base_case_trial = Trial(base_case_params)
base_case_trial.run_trial()
base_case_trial.calculate_trial_results()
what_if_params = Param(num_nurses=2, num_receptionists=2)
what_if_trial = Trial(what_if_params)
what_if_trial.run_trial()
what_if_trial.calculate_trial_results()
print ("BASE CASE TRIAL RESULTS")
print ("-----------------------")
print ("Queuing Time for Registration")
print (f"Mean : {base_case_trial.trial_mean_q_time_registration:.2f} minutes")
print (f"SD : {base_case_trial.trial_sd_q_time_registration:.2f} minutes")
print (
f"90th Perc : {base_case_trial.trial_perc_90_q_time_registration:.2f}",
"minutes"
)
print (f"Standard Error : {base_case_trial.se_q_time_registration:.2f}")
print (
f"95% CI : ({base_case_trial.ci_lower_q_time_registration:.2f}, ",
f"{base_case_trial.ci_upper_q_time_registration:.2f}) minutes"
)
print ()
print ("Queuing Time for the Nurse")
print (f"Mean : {base_case_trial.trial_mean_q_time_nurse:.2f} minutes")
print (f"SD : {base_case_trial.trial_sd_q_time_nurse:.2f} minutes")
print (f"90th Perc : {base_case_trial.trial_perc_90_q_time_nurse:.2f} minutes")
print (f"Standard Error : {base_case_trial.se_q_time_nurse:.2f}")
print (
f"95% CI : ({base_case_trial.ci_lower_q_time_nurse:.2f}, ",
f"{base_case_trial.ci_upper_q_time_nurse:.2f}) minutes"
)
print ()
# WRITE YOUR NEW CODE BELOW
print ("2 NURSES, 2 RECEPTIONISTS TRIAL RESULTS")
print ("----------------------")
print ("Queuing Time for Registration")
print (f"Mean : {what_if_trial.trial_mean_q_time_registration:.2f} minutes")
print (f"SD : {what_if_trial.trial_sd_q_time_registration:.2f} minutes")
print (
f"90th Perc : {what_if_trial.trial_perc_90_q_time_registration:.2f}",
"minutes"
)
print (f"Standard Error : {what_if_trial.se_q_time_registration:.2f}")
print (
f"95% CI : ({what_if_trial.ci_lower_q_time_registration:.2f}, ",
f"{what_if_trial.ci_upper_q_time_registration:.2f}) minutes"
)
print ()
print ("Queuing Time for the Nurse")
print (f"Mean : {what_if_trial.trial_mean_q_time_nurse:.2f} minutes")
print (f"SD : {what_if_trial.trial_sd_q_time_nurse:.2f} minutes")
print (f"90th Perc : {what_if_trial.trial_perc_90_q_time_nurse:.2f} minutes")
print (f"Standard Error : {what_if_trial.se_q_time_nurse:.2f}")
print (
f"95% CI : ({what_if_trial.ci_lower_q_time_nurse:.2f}, ",
f"{what_if_trial.ci_upper_q_time_nurse:.2f}) minutes"
)
print ()
# WRITE YOUR NEW CODE BELOW