Advanced trace#
1. Imports#
The model used to demonstrate trace functionality is implemented in simpy
import simpy
simpy.__version__
'4.1.1'
import sim_tools
sim_tools.__version__
'0.7.1'
from sim_tools.distributions import (
Exponential,
Normal,
Lognormal,
Bernoulli
)
from sim_tools.time_dependent import NSPPThinning
from sim_tools.trace import Traceable
import numpy as np
import pandas as pd
import itertools
2. Constants and defaults for modelling as-is#
2.1 Distribution parameters#
# sign-in/triage parameters
DEFAULT_TRIAGE_MEAN = 3.0
# registration parameters
DEFAULT_REG_MEAN = 5.0
DEFAULT_REG_VAR= 2.0
# examination parameters
DEFAULT_EXAM_MEAN = 16.0
DEFAULT_EXAM_VAR = 3.0
DEFAULT_EXAM_MIN = 0.5
# trauma/stabilisation
DEFAULT_TRAUMA_MEAN = 90.0
# Trauma treatment
DEFAULT_TRAUMA_TREAT_MEAN = 30.0
DEFAULT_TRAUMA_TREAT_VAR = 4.0
# Non trauma treatment
DEFAULT_NON_TRAUMA_TREAT_MEAN = 13.3
DEFAULT_NON_TRAUMA_TREAT_VAR = 2.0
# prob patient requires treatment given trauma
DEFAULT_NON_TRAUMA_TREAT_P = 0.60
# proportion of patients triaged as trauma
DEFAULT_PROB_TRAUMA = 0.12
2.2 Time dependent arrival rates data#
The data for arrival rates varies between clinic opening at 6am and closure at 12am.
# data are held in the Github repo and loaded from there.
NSPP_PATH = 'https://raw.githubusercontent.com/TomMonks/' \
+ 'open-science-for-sim/main/src/notebooks/01_foss_sim/data/ed_arrivals.csv'
2.3 Resource counts#
Inter count variables representing the number of resources at each activity in the processes.
DEFAULT_N_TRIAGE = 1
DEFAULT_N_REG = 1
DEFAULT_N_EXAM = 3
DEFAULT_N_TRAUMA = 2
# Non-trauma cubicles
DEFAULT_N_CUBICLES_1 = 1
# trauma pathway cubicles
DEFAULT_N_CUBICLES_2 = 1
2.4 Simulation model run settings#
# default random number SET
DEFAULT_RNG_SET = None
N_STREAMS = 20
# default results collection period
DEFAULT_RESULTS_COLLECTION_PERIOD = 60 * 19
# number of replications.
DEFAULT_N_REPS = 5
## set to True to show a trace of the simulation model (best used in single run mode)
DEBUG = False
# provide a list of integers that represent patient IDs. Only these patients will be tracked
TRACKED = None
3. Model parameterisation#
For convienience a container class is used to hold the large number of model parameters. The Scenario class includes defaults these can easily be changed and at runtime to experiments with different designs.
def load_arrival_profile(profile_path):
df = (pd.read_csv(profile_path)
.reset_index()
.rename(columns={'index':'t'})
.assign(mean_iat = lambda x: 60.0 / x.arrival_rate,
t = lambda x: x.t * 60.0)
)
return df
load_arrival_profile(NSPP_PATH)
| t | period | arrival_rate | mean_iat | |
|---|---|---|---|---|
| 0 | 0.0 | 6AM-7AM | 2.366667 | 25.352113 |
| 1 | 60.0 | 7AM-8AM | 2.800000 | 21.428571 |
| 2 | 120.0 | 8AM-9AM | 8.833333 | 6.792453 |
| 3 | 180.0 | 9AM-10AM | 10.433333 | 5.750799 |
| 4 | 240.0 | 10AM-11AM | 14.800000 | 4.054054 |
| 5 | 300.0 | 11AM-12PM | 26.266667 | 2.284264 |
| 6 | 360.0 | 12PM-1PM | 31.400000 | 1.910828 |
| 7 | 420.0 | 1PM-2PM | 18.066667 | 3.321033 |
| 8 | 480.0 | 2PM-3PM | 16.466667 | 3.643725 |
| 9 | 540.0 | 3PM-4PM | 12.033333 | 4.986150 |
| 10 | 600.0 | 4PM-5PM | 11.600000 | 5.172414 |
| 11 | 660.0 | 5PM-6PM | 28.866667 | 2.078522 |
| 12 | 720.0 | 6PM-7PM | 18.033333 | 3.327172 |
| 13 | 780.0 | 7PM-8PM | 11.500000 | 5.217391 |
| 14 | 840.0 | 8PM-9PM | 5.300000 | 11.320755 |
| 15 | 900.0 | 9PM-10PM | 4.066667 | 14.754098 |
| 16 | 960.0 | 10PM-11PM | 2.200000 | 27.272727 |
| 17 | 1020.0 | 11PM-12AM | 2.100000 | 28.571429 |
class Scenario:
'''
Container class for scenario parameters/arguments
Passed to a model and its process classes
'''
def __init__(self, random_number_set=DEFAULT_RNG_SET,
n_triage=DEFAULT_N_TRIAGE,
n_reg=DEFAULT_N_REG,
n_exam=DEFAULT_N_EXAM,
n_trauma=DEFAULT_N_TRAUMA,
n_cubicles_1=DEFAULT_N_CUBICLES_1,
n_cubicles_2=DEFAULT_N_CUBICLES_2,
triage_mean=DEFAULT_TRIAGE_MEAN,
reg_mean=DEFAULT_REG_MEAN,
reg_var=DEFAULT_REG_VAR,
exam_mean=DEFAULT_EXAM_MEAN,
exam_var=DEFAULT_EXAM_VAR,
exam_min=DEFAULT_EXAM_MIN,
trauma_mean=DEFAULT_TRAUMA_MEAN,
trauma_treat_mean=DEFAULT_TRAUMA_TREAT_MEAN,
trauma_treat_var=DEFAULT_TRAUMA_TREAT_VAR,
non_trauma_treat_mean=DEFAULT_NON_TRAUMA_TREAT_MEAN,
non_trauma_treat_var=DEFAULT_NON_TRAUMA_TREAT_VAR,
non_trauma_treat_p=DEFAULT_NON_TRAUMA_TREAT_P,
prob_trauma=DEFAULT_PROB_TRAUMA,
debug=DEBUG,
tracked=TRACKED):
'''
Create a scenario to parameterise the simulation model
Parameters:
-----------
random_number_set: int, optional (default=DEFAULT_RNG_SET)
Set to control the initial seeds of each stream of pseudo
random numbers used in the model.
n_triage: int
The number of triage cubicles
n_reg: int
The number of registration clerks
n_exam: int
The number of examination rooms
n_trauma: int
The number of trauma bays for stabilisation
n_cubicles_1: int
The number of non-trauma treatment cubicles
n_cubicles_2: int
The number of trauma treatment cubicles
triage_mean: float
Mean duration of the triage distribution (Exponential)
reg_mean: float
Mean duration of the registration distribution (Lognormal)
reg_var: float
Variance of the registration distribution (Lognormal)
exam_mean: float
Mean of the examination distribution (Normal)
exam_var: float
Variance of the examination distribution (Normal)
exam_min: float
The minimum value that an examination can take (Truncated Normal)
trauma_mean: float
Mean of the trauma stabilisation distribution (Exponential)
trauma_treat_mean: float
Mean of the trauma cubicle treatment distribution (Lognormal)
trauma_treat_var: float
Variance of the trauma cubicle treatment distribution (Lognormal)
non_trauma_treat_mean: float
Mean of the non trauma treatment distribution
non_trauma_treat_var: float
Variance of the non trauma treatment distribution
non_trauma_treat_p: float
Probability non trauma patient requires treatment
prob_trauma: float
probability that a new arrival is a trauma patient.
debug: int
Set to true to display a trace of simulated events in the model
Note this is best used in single_run mode.
tracked: list
List of integers that represent patient IDs.
Used with debug. If debug is True then only display trace for specific
patients.
'''
# sampling
self.random_number_set = random_number_set
# store parameters for sampling
self.triage_mean = triage_mean
self.reg_mean = reg_mean
self.reg_var = reg_var
self.exam_mean= exam_mean
self.exam_var = exam_var
self.exam_min = exam_min
self.trauma_mean = trauma_mean
self.trauma_treat_mean = trauma_treat_mean
self.trauma_treat_var = trauma_treat_var
self.non_trauma_treat_mean = non_trauma_treat_mean
self.non_trauma_treat_var = non_trauma_treat_var
self.non_trauma_treat_p = non_trauma_treat_p
self.prob_trauma = prob_trauma
self.init_sampling()
# count of each type of resource
self.init_resourse_counts(n_triage, n_reg, n_exam, n_trauma,
n_cubicles_1, n_cubicles_2)
# debug/trace information
self.debug = debug
self.tracked = tracked
def set_random_no_set(self, random_number_set):
'''
Controls the random sampling
Parameters:
----------
random_number_set: int
Used to control the set of pseudo random numbers
used by the distributions in the simulation.
'''
self.random_number_set = random_number_set
self.init_sampling()
def init_resourse_counts(self, n_triage, n_reg, n_exam, n_trauma,
n_cubicles_1, n_cubicles_2):
'''
Init the counts of resources to default values...
'''
self.n_triage = n_triage
self.n_reg = n_reg
self.n_exam = n_exam
self.n_trauma = n_trauma
# non-trauma (1), trauma (2) treatment cubicles
self.n_cubicles_1 = n_cubicles_1
self.n_cubicles_2 = n_cubicles_2
def init_sampling(self):
'''
Create the distributions used by the model and initialise
the random seeds of each.
'''
# MODIFICATION. Better method for producing n non-overlapping streams
seed_sequence = np.random.SeedSequence(self.random_number_set)
# Generate n high quality child seeds
self.seeds = seed_sequence.spawn(N_STREAMS)
# create distributions
# Triage duration
self.triage_dist = Exponential(self.triage_mean,
random_seed=self.seeds[0])
# Registration duration (non-trauma only)
self.reg_dist = Lognormal(self.reg_mean,
np.sqrt(self.reg_var),
random_seed=self.seeds[1])
# Evaluation (non-trauma only)
self.exam_dist = Normal(self.exam_mean,
np.sqrt(self.exam_var),
minimum=self.exam_min,
random_seed=self.seeds[2])
# Trauma/stablisation duration (trauma only)
self.trauma_dist = Exponential(self.trauma_mean,
random_seed=self.seeds[3])
# Non-trauma treatment
self.nt_treat_dist = Lognormal(self.non_trauma_treat_mean,
np.sqrt(self.non_trauma_treat_var),
random_seed=self.seeds[4])
# treatment of trauma patients
self.treat_dist = Lognormal(self.trauma_treat_mean,
np.sqrt(self.trauma_treat_var),
random_seed=self.seeds[5])
# probability of non-trauma patient requiring treatment
self.nt_p_treat_dist = Bernoulli(self.non_trauma_treat_p,
random_seed=self.seeds[6])
# probability of non-trauma versus trauma patient
self.p_trauma_dist = Bernoulli(self.prob_trauma,
random_seed=self.seeds[7])
# init sampling for non-stationary poisson process
self.init_nspp()
def init_nspp(self):
# read arrival profile
self.arrivals = load_arrival_profile(NSPP_PATH)
# maximum arrival rate (smallest time between arrivals)
self.lambda_max = self.arrivals['arrival_rate'].max()
# NSPP thinning distribution
self.thinning_dist = NSPPThinning(self.arrivals,
self.seeds[8],
self.seeds[9])
6. Patient Pathways Process Logic#
simpy uses a process based worldview. We can easily create whatever logic - simple or complex for the model. Here the process logic for trauma and non-trauma patients is seperated into two classes TraumaPathway and NonTraumaPathway.
class TraumaPathway(Traceable):
'''
Encapsulates the process a patient with severe injuries or illness.
These patients are signed into the ED and triaged as having severe injuries
or illness.
Patients are stabilised in resus (trauma) and then sent to Treatment.
Following treatment they are discharged.
'''
def __init__(self, identifier, env, args):
'''
Constructor method
Params:
-----
identifier: int
a numeric identifier for the patient.
env: simpy.Environment
the simulation environment
args: Scenario
Container class for the simulation parameters
'''
# initialise Trace
super().__init__(debug=args.debug)
self.identifier = identifier
self.env = env
self.args = args
# metrics
self.arrival = -np.inf
self.wait_triage = -np.inf
self.wait_trauma = -np.inf
self.wait_treat = -np.inf
self.total_time = -np.inf
self.triage_duration = -np.inf
self.trauma_duration = -np.inf
self.treat_duration = -np.inf
def execute(self):
'''
simulates the major treatment process for a patient
1. request and wait for sign-in/triage
2. trauma
3. treatment
'''
# record the time of arrival and entered the triage queue
self.arrival = self.env.now
self.trauma_arrival()
# request sign-in/triage
with self.args.triage.request() as req:
yield req
# record the waiting time for triage
self.wait_triage = self.env.now - self.arrival
# triage begin event
self.triage_begin()
# sample triage duration.
self.triage_duration = self.args.triage_dist.sample()
yield self.env.timeout(self.triage_duration)
# triage complete event
self.triage_complete()
# record the time that entered the trauma queue
start_wait = self.env.now
# request trauma room
with self.args.trauma.request() as req:
yield req
# record the waiting time for trauma
self.wait_trauma = self.env.now - start_wait
# trauma begin event
self.trauma_begin()
# sample stablisation duration.
self.trauma_duration = self.args.trauma_dist.sample()
yield self.env.timeout(self.trauma_duration)
# trauma complete event
self.trauma_complete()
# record the time that entered the treatment queue
start_wait = self.env.now
# request treatment cubicle
with self.args.cubicle_2.request() as req:
yield req
# record the waiting time for trauma
self.wait_treat = self.env.now - start_wait
self.treatment_begin()
# sample treatment duration.
self.treat_duration = self.args.treat_dist.sample()
yield self.env.timeout(self.treat_duration)
self.treatment_complete()
# total time in system
self.total_time = self.env.now - self.arrival
def trauma_arrival(self):
'''Trauma arrival event'''
self.trace(self.env.now, 'arrival at centre 🚑', self.identifier)
def triage_begin(self):
'''Enter triage event'''
self.trace(self.env.now, f'enter triage. Waiting time: {self.wait_triage:.3f}',
self.identifier)
def trauma_begin(self):
'''Enter trauma stabilisation'''
self.trace(self.env.now, f'enter stabilisation. Waiting time: {self.wait_trauma:.3f}',
self.identifier)
def treatment_begin(self):
'''Enter trauma treatment post stabilisation'''
self.trace(self.env.now, f'enter treatment. Waiting time: {self.wait_treat:.3f}',
self.identifier)
def triage_complete(self):
'''Triage complete event'''
self.trace(self.env.now, 'triage complete', self.identifier)
def trauma_complete(self):
'''Patient stay in trauma is complete.'''
self.trace(self.env.now, 'stabilisation complete.', self.identifier)
def treatment_complete(self):
'''Treatment complete event'''
self.trace(self.env.now, f'patient {self.identifier} treatment complete; '
f'waiting time was {self.wait_treat:.3f}', self.identifier)
def _trace_config(self):
'''Override trace config'''
config = {
"name":"Trauma",
"name_colour":"bold magenta",
"time_colour":'bold magenta',
"tracked":self.args.tracked
}
return config
class NonTraumaPathway(Traceable):
'''
Encapsulates the process a patient with minor injuries and illness.
These patients are signed into the ED and triaged as having minor
complaints and streamed to registration and then examination.
Post examination 40% are discharged while 60% proceed to treatment.
Following treatment they are discharged.
'''
def __init__(self, identifier, env, args):
'''
Constructor method
Params:
-----
identifier: int
a numeric identifier for the patient.
env: simpy.Environment
the simulation environment
args: Scenario
Container class for the simulation parameters
'''
super().__init__(debug=args.debug)
self.identifier = identifier
self.env = env
self.args = args
# triage resource
self.triage = args.triage
# metrics
self.arrival = -np.inf
self.wait_triage = -np.inf
self.wait_reg = -np.inf
self.wait_exam = -np.inf
self.wait_treat = -np.inf
self.total_time = -np.inf
self.triage_duration = -np.inf
self.reg_duration = -np.inf
self.exam_duration = -np.inf
self.treat_duration = -np.inf
def execute(self):
'''
simulates the non-trauma/minor treatment process for a patient
1. request and wait for sign-in/triage
2. patient registration
3. examination
4. split patients
4.1 % discharged
4.2 % treatment then discharge
'''
# record the time of arrival and entered the triage queue
self.arrival = self.env.now
# trace arrival
self.trace(self.env.now, 'arrival to centre.', self.identifier)
# request sign-in/triage
with self.triage.request() as req:
yield req
# record the waiting time for triage
self.wait_triage = self.env.now - self.arrival
# trace enter triage
self.trace(self.env.now, 'entering triage. '
f'Waiting time: {self.wait_triage:.3f} mins',
self.identifier)
# sample triage duration.
self.triage_duration = self.args.triage_dist.sample()
yield self.env.timeout(self.triage_duration)
# trace exit triage
self.trace(self.env.now, 'triage complete', self.identifier)
# record the time that entered the registration queue
start_wait = self.env.now
# request registration clert
with self.args.registration.request() as req:
yield req
# record the waiting time for registration
self.wait_reg = self.env.now - start_wait
# trace begin registration
self.trace(self.env.now, 'starting patient registration.', self.identifier)
# sample registration duration.
self.reg_duration = self.args.reg_dist.sample()
yield self.env.timeout(self.reg_duration)
# registration complete...
self.trace(self.env.now, 'patient registered;'
f'waiting time was {self.wait_triage:.3f}',
self.identifier)
# record the time that entered the evaluation queue
start_wait = self.env.now
# request examination resource
with self.args.exam.request() as req:
yield req
# record the waiting time for registration
self.wait_exam = self.env.now - start_wait
# trace begin examination
self.trace(self.env.now, f'enter examination. Waiting time: {self.wait_exam:.3f}',
self.identifier)
# sample examination duration.
self.exam_duration = self.args.exam_dist.sample()
yield self.env.timeout(self.exam_duration)
# trace exit examination
self.trace(self.env.now, 'examination complete.', self.identifier)
# sample if patient requires treatment?
self.require_treat = self.args.nt_p_treat_dist.sample()
if self.require_treat:
# record the time that entered the treatment queue
start_wait = self.env.now
# request treatment cubicle
with self.args.cubicle_1.request() as req:
yield req
# record the waiting time for treatment
self.wait_treat = self.env.now - start_wait
# trace enter treatment
self.trace(self.env.now, f'enter treatment. Waiting time:{self.wait_treat:.3f}',
self.identifier)
# sample treatment duration.
self.treat_duration = self.args.nt_treat_dist.sample()
yield self.env.timeout(self.treat_duration)
# trace exit treatment
self.trace(self.env.now, 'treatment complete â›”', self.identifier)
# total time in system
self.total_time = self.env.now - self.arrival
def _trace_config(self):
'''Override trace config'''
config = {
"name":"Non-trauma",
"name_colour":"bold green",
"time_colour":'bold green',
"message_colour":'black',
"tracked":self.args.tracked
}
return config
7. Main model class#
The main class that a user interacts with to run the model is TreatmentCentreModel. This implements a .run() method, contains a simple algorithm for the non-stationary poission process for patients arrivals and inits instances of TraumaPathway or NonTraumaPathway depending on the arrival type.
class TreatmentCentreModel:
'''
The treatment centre model
Patients arrive at random to a treatment centre, are triaged
and then processed in either a trauma or non-trauma pathway.
'''
def __init__(self, args):
self.env = simpy.Environment()
self.args = args
self.init_resources()
self.patients = []
self.trauma_patients = []
self.non_trauma_patients = []
self.rc_period = None
self.results = None
def init_resources(self):
'''
Init the number of resources
and store in the arguments container object
Resource list:
1. Sign-in/triage bays
2. registration clerks
3. examination bays
4. trauma bays
5. non-trauma cubicles (1)
6. trauma cubicles (2)
'''
# sign/in triage
self.args.triage = simpy.Resource(self.env,
capacity=self.args.n_triage)
# registration
self.args.registration = simpy.Resource(self.env,
capacity=self.args.n_reg)
# examination
self.args.exam = simpy.Resource(self.env,
capacity=self.args.n_exam)
# trauma
self.args.trauma = simpy.Resource(self.env,
capacity=self.args.n_trauma)
# non-trauma treatment
self.args.cubicle_1 = simpy.Resource(self.env,
capacity=self.args.n_cubicles_1)
# trauma treatment
self.args.cubicle_2 = simpy.Resource(self.env,
capacity=self.args.n_cubicles_2)
def run(self, results_collection_period=DEFAULT_RESULTS_COLLECTION_PERIOD):
'''
Conduct a single run of the model in its current
configuration
Parameters:
----------
results_collection_period, float, optional
default = DEFAULT_RESULTS_COLLECTION_PERIOD
warm_up, float, optional (default=0)
length of initial transient period to truncate
from results.
Returns:
--------
None
'''
# setup the arrival generator process
self.env.process(self.arrivals_generator())
# store rc perio
self.rc_period = results_collection_period
# run
self.env.run(until=results_collection_period)
def arrivals_generator(self):
'''
Simulate the arrival of patients to the model
Patients either follow a TraumaPathway or
NonTraumaPathway simpy process.
Non stationary arrivals implemented via Thinning acceptance-rejection
algorithm.
'''
for patient_count in itertools.count():
# sample iat via thinning
interarrival_time = self.args.thinning_dist.sample(self.env.now)
# delay until next arrival
yield self.env.timeout(interarrival_time)
# sample if the patient is trauma or non-trauma
trauma = self.args.p_trauma_dist.sample()
if trauma:
# create and store a trauma patient to update KPIs.
new_patient = TraumaPathway(patient_count, self.env, self.args)
self.trauma_patients.append(new_patient)
else:
# create and store a non-trauma patient to update KPIs.
new_patient = NonTraumaPathway(patient_count, self.env, self.args)
self.non_trauma_patients.append(new_patient)
# start the pathway process for the patient
self.env.process(new_patient.execute())
8. Logic to process end of run results.#
the class SimulationSummary accepts a TraumaCentreModel. At the end of a run it can be used calculate mean queuing times and the percentage of the total run that a resource was in use.
class SimulationSummary:
'''
End of run result processing logic of the simulation model
'''
def __init__(self, model):
'''
Constructor
Params:
------
model: TraumaCentreModel
The model.
'''
self.model = model
self.args = model.args
self.results = None
def process_run_results(self):
'''
Calculates statistics at end of run.
'''
self.results = {}
# list of all patients
patients = self.model.non_trauma_patients + self.model.trauma_patients
# mean triage times (both types of patient)
mean_triage_wait = self.get_mean_metric('wait_triage', patients)
# triage utilisation (both types of patient)
triage_util = self.get_resource_util('triage_duration',
self.args.n_triage,
patients)
# mean waiting time for registration (non_trauma)
mean_reg_wait = self.get_mean_metric('wait_reg',
self.model.non_trauma_patients)
# registration utilisation (trauma)
reg_util = self.get_resource_util('reg_duration',
self.args.n_reg,
self.model.non_trauma_patients)
# mean waiting time for examination (non_trauma)
mean_wait_exam = self.get_mean_metric('wait_exam',
self.model.non_trauma_patients)
# examination utilisation (non-trauma)
exam_util = self.get_resource_util('exam_duration',
self.args.n_exam,
self.model.non_trauma_patients)
# mean waiting time for treatment (non-trauma)
mean_treat_wait = self.get_mean_metric('wait_treat',
self.model.non_trauma_patients)
# treatment utilisation (non_trauma)
treat_util1 = self.get_resource_util('treat_duration',
self.args.n_cubicles_1,
self.model.non_trauma_patients)
# mean total time (non_trauma)
mean_total = self.get_mean_metric('total_time',
self.model.non_trauma_patients)
# mean waiting time for trauma
mean_trauma_wait = self.get_mean_metric('wait_trauma',
self.model.trauma_patients)
# trauma utilisation (trauma)
trauma_util = self.get_resource_util('trauma_duration',
self.args.n_trauma,
self.model.trauma_patients)
# mean waiting time for treatment (rauma)
mean_treat_wait2 = self.get_mean_metric('wait_treat',
self.model.trauma_patients)
# treatment utilisation (trauma)
treat_util2 = self.get_resource_util('treat_duration',
self.args.n_cubicles_2,
self.model.trauma_patients)
# mean total time (trauma)
mean_total2 = self.get_mean_metric('total_time',
self.model.trauma_patients)
self.results = {'00_arrivals':len(patients),
'01a_triage_wait': mean_triage_wait,
'01b_triage_util': triage_util,
'02a_registration_wait':mean_reg_wait,
'02b_registration_util': reg_util,
'03a_examination_wait':mean_wait_exam,
'03b_examination_util': exam_util,
'04a_treatment_wait(non_trauma)':mean_treat_wait,
'04b_treatment_util(non_trauma)':treat_util1,
'05_total_time(non-trauma)':mean_total,
'06a_trauma_wait':mean_trauma_wait,
'06b_trauma_util':trauma_util,
'07a_treatment_wait(trauma)':mean_treat_wait2,
'07b_treatment_util(trauma)':treat_util2,
'08_total_time(trauma)':mean_total2,
'09_throughput': self.get_throughput(patients)}
def get_mean_metric(self, metric, patients):
'''
Calculate mean of the performance measure for the
select cohort of patients,
Only calculates metrics for patients where it has been
measured.
Params:
-------
metric: str
The name of the metric e.g. 'wait_treat'
patients: list
A list of patients
'''
mean = np.array([getattr(p, metric) for p in patients
if getattr(p, metric) > -np.inf]).mean()
return mean
def get_resource_util(self, metric, n_resources, patients):
'''
Calculate proportion of the results collection period
where a resource was in use.
Done by tracking the duration by patient.
Only calculates metrics for patients where it has been
measured.
Params:
-------
metric: str
The name of the metric e.g. 'treatment_duration'
patients: list
A list of patients
'''
total = np.array([getattr(p, metric) for p in patients
if getattr(p, metric) > -np.inf]).sum()
return total / (self.model.rc_period * n_resources)
def get_throughput(self, patients):
'''
Returns the total number of patients that have successfully
been processed and discharged in the treatment centre
(they have a total time record)
Params:
-------
patients: list
list of all patient objects simulated.
Returns:
------
float
'''
return len([p for p in patients if p.total_time > -np.inf])
def summary_frame(self):
'''
Returns run results as a pandas.DataFrame
Returns:
-------
pd.DataFrame
'''
#append to results df
if self.results is None:
self.process_run_results()
df = pd.DataFrame({'1':self.results})
df = df.T
df.index.name = 'rep'
return df
9. Model execution#
We note that there are many ways to setup a simpy model and execute it (that is part of its fantastic flexibility). The organisation of code we show below is based on our experience of using the package in practice. The approach also allows for easy parallisation over multiple CPU cores using joblib.
We include two functions. single_run() and multiple_replications. The latter is used to repeatedly call and process the results from single_run.
def single_run(scenario, rc_period=DEFAULT_RESULTS_COLLECTION_PERIOD,
random_no_set=DEFAULT_RNG_SET):
'''
Perform a single run of the model and return the results
Parameters:
-----------
scenario: Scenario object
The scenario/paramaters to run
rc_period: int
The length of the simulation run that collects results
random_no_set: int or None, optional (default=DEFAULT_RNG_SET)
Controls the set of random seeds used by the stochastic parts of the
model. Set to different ints to get different results. Set to None
for a random set of seeds.
Returns:
--------
pandas.DataFrame:
results from single run.
'''
# set random number set - this controls sampling for the run.
scenario.set_random_no_set(random_no_set)
# create an instance of the model
model = TreatmentCentreModel(scenario)
# run the model
model.run(results_collection_period=rc_period)
# run results
summary = SimulationSummary(model)
summary_df = summary.summary_frame()
return summary_df
def multiple_replications(scenario, rc_period=DEFAULT_RESULTS_COLLECTION_PERIOD,
n_reps=5):
'''
Perform multiple replications of the model.
Params:
------
scenario: Scenario
Parameters/arguments to configurethe model
rc_period: float, optional (default=DEFAULT_RESULTS_COLLECTION_PERIOD)
results collection period.
the number of minutes to run the model to collect results
n_reps: int, optional (default=DEFAULT_N_REPS)
Number of independent replications to run.
Returns:
--------
pandas.DataFrame
'''
results = [single_run(scenario, rc_period, random_no_set=rep)
for rep in range(n_reps)]
#format and return results in a dataframe
df_results = pd.concat(results)
df_results.index = np.arange(1, len(df_results)+1)
df_results.index.name = 'rep'
return df_results
9.1 Single run of the model#
The script below performs a single replication of the simulation model.
Try:
Changing the
random_no_setof thesingle_runcall.Assigning the value
TruetodebugTracking patient arrivals 10 and 19 through the model by setting
tracked=[10, 19]
# create the default scenario.
args = Scenario(debug=True, tracked=[10, 19])
# use the single_run() func
# try changing `random_no_set` to see different run results
print('Running simulation ...', end=' => ')
results = single_run(args, random_no_set=42)
print('simulation complete.')
# show results (transpose replication for easier view)
results.T
Running simulation ... =>
[169.77]:arrival to centre.
[175.62]:entering triage. Waiting time: 5.849 mins
[175.79]:triage complete
[188.70]:starting patient registration.
[193.51]:patient registered;waiting time was 5.849
[193.51]:enter examination. Waiting time: 0.000
[204.96]:arrival at centre 🚑
[208.94]:examination complete.
[208.94]:enter treatment. Waiting time:0.000
[214.79]:enter triage. Waiting time: 9.833
[221.03]:triage complete
[221.63]:treatment complete â›”
[384.82]:enter stabilisation. Waiting time: 163.799
[441.65]:stabilisation complete.
[473.10]:enter treatment. Waiting time: 31.447
[505.03]:patient 19 treatment complete; waiting time was 31.447
simulation complete.
| rep | 1 |
|---|---|
| 00_arrivals | 209.000000 |
| 01a_triage_wait | 16.622674 |
| 01b_triage_util | 0.527512 |
| 02a_registration_wait | 111.161345 |
| 02b_registration_util | 0.801061 |
| 03a_examination_wait | 24.927965 |
| 03b_examination_util | 0.851285 |
| 04a_treatment_wait(non_trauma) | 172.435861 |
| 04b_treatment_util(non_trauma) | 0.845652 |
| 05_total_time(non-trauma) | 248.848441 |
| 06a_trauma_wait | 201.144403 |
| 06b_trauma_util | 0.919830 |
| 07a_treatment_wait(trauma) | 22.620880 |
| 07b_treatment_util(trauma) | 0.493576 |
| 08_total_time(trauma) | 310.384648 |
| 09_throughput | 155.000000 |
9.2 Disabling trace for multiple independent replications#
Important: set debug=False in Scenario for multiple replications. This will ensure that each replication is not traced.
args = Scenario(debug=False)
#run multiple replications.
results = multiple_replications(args, n_reps=50)
print("All replications complete.")
# summarise the results (2.dp)
results.mean().round(2)
All replications complete.
00_arrivals 227.72
01a_triage_wait 35.24
01b_triage_util 0.61
02a_registration_wait 105.57
02b_registration_util 0.84
03a_examination_wait 25.55
03b_examination_util 0.85
04a_treatment_wait(non_trauma) 136.66
04b_treatment_util(non_trauma) 0.87
05_total_time(non-trauma) 234.34
06a_trauma_wait 151.68
06b_trauma_util 0.83
07a_treatment_wait(trauma) 14.31
07b_treatment_util(trauma) 0.50
08_total_time(trauma) 292.28
09_throughput 162.16
dtype: float64