import io
import numpy as np
import pandas as pd
import plotly.io as pio
import simpy
from sim_tools.distributions import Exponential, Lognormal
from vidigi.logging import EventLogger, TrialLogger
from vidigi.resources import VidigiStore
from vidigi.utils import EventPosition, create_event_position_df
pio.renderers.default = "notebook"Feature Example: The TrialLogger Class
While we could put the list of EventLoggers from each run of a trial into a plain list, wrapping them in a TrialLogger instead unlocks a set of helper methods for statistics and plots computed across every run at once - which is most of what’s new in vidigi 2.0.0.
This notebook picks up where the EventLogger notebook leaves off, so see that one first for a walkthrough of the EventLogger class itself, event position dataframes, and animation. We reuse the same multi-resource, branching clinic model here (collapsed below, since it’s already explained there), since it gives TrialLogger’s plots more interesting data to work with than the simpler single-resource example.
Model setup
Show the import code
# Import additional required distributions
from sim_tools.distributions import Bernoulli, Normal, UniformShow the global parameter class code
# 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:
"""
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 stablisation
n_cubicles_non_trauma_treat: int
The number of non-trauma treatment cubicles
n_cubicles_trauma_treat: 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)
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.
"""
random_number_set = 42
n_triage = 2
n_reg = 2
n_exam = 3
n_trauma = 4
n_cubicles_non_trauma_treat = 4
n_cubicles_trauma_treat = 5
triage_mean = 6
reg_mean = 8
reg_var = 2
exam_mean = 16
exam_var = 3
trauma_mean = 90
trauma_treat_mean = 30
trauma_treat_var = 4
non_trauma_treat_mean = 13.3
non_trauma_treat_var = 2
non_trauma_treat_p = 0.6
prob_trauma = 0.12
arrival_df = "ed_arrivals.csv"
sim_duration = 600
number_of_runs = 100Show the patient class code
class Patient:
"""
Class defining details for a patient entity
"""
def __init__(self, p_id):
"""
Constructor method
Params:
-----
identifier: int
a numeric identifier for the patient.
"""
self.identifier = p_id
# Time of arrival in model/at centre
self.arrival = -np.inf
# Total time in pathway
self.total_time = -np.inf
# Shared waits
self.wait_triage = -np.inf
self.wait_reg = -np.inf
self.wait_treat = -np.inf
# Non-trauma pathway - examination wait
self.wait_exam = -np.inf
# Trauma pathway - stabilisation wait
self.wait_trauma = -np.inf
# Shared durations
self.triage_duration = -np.inf
self.reg_duration = -np.inf
self.treat_duration = -np.inf
# Non-trauma pathway - examination duration
self.exam_duration = -np.inf
# Trauma pathway - stabilisation duration
self.trauma_duration = -np.infShow the model code
# Class representing our model of the clinic.
class Model:
"""
Simulates the simplest minor treatment process for a patient
1. Arrive
2. Examined/treated by nurse when one available
3. Discharged
"""
# 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, n_cubicles_trauma_treat=None):
# Create a SimPy environment in which everything will live
self.env = simpy.Environment()
# Store the passed in run number
self.run_number = run_number
self.logger = EventLogger(env=self.env, run_number=self.run_number)
# n_cubicles_trauma_treat= is an override for the scenario-comparison
# section further down; every other cell in this notebook omits it, so
# self.n_cubicles_trauma_treat is a verified no-op equal to
# g.n_cubicles_trauma_treat.
self.n_cubicles_trauma_treat = (
n_cubicles_trauma_treat
if n_cubicles_trauma_treat is not None
else g.n_cubicles_trauma_treat
)
# Create a patient counter (which we'll use as a patient ID)
self.patient_counter = 0
self.trauma_patients = []
self.non_trauma_patients = []
# Create our resources
self.init_resources()
# Create our distributions
self.init_distributions()
def init_distributions(self):
# Create distributions
# Triage duration
self.triage_dist = Exponential(
g.triage_mean, random_seed=self.run_number * g.random_number_set
)
# Registration duration (non-trauma only)
self.reg_dist = Lognormal(
g.reg_mean,
np.sqrt(g.reg_var),
random_seed=self.run_number * g.random_number_set,
)
# Evaluation (non-trauma only)
self.exam_dist = Normal(
g.exam_mean,
np.sqrt(g.exam_var),
random_seed=self.run_number * g.random_number_set,
)
# Trauma/stablisation duration (trauma only)
self.trauma_dist = Exponential(
g.trauma_mean, random_seed=self.run_number * g.random_number_set
)
# Non-trauma treatment
self.nt_treat_dist = Lognormal(
g.non_trauma_treat_mean,
np.sqrt(g.non_trauma_treat_var),
random_seed=self.run_number * g.random_number_set,
)
# treatment of trauma patients
self.treat_dist = Lognormal(
g.trauma_treat_mean,
np.sqrt(g.non_trauma_treat_var),
random_seed=self.run_number * g.random_number_set,
)
# probability of non-trauma patient requiring treatment
self.nt_p_treat_dist = Bernoulli(
g.non_trauma_treat_p, random_seed=self.run_number * g.random_number_set
)
# probability of non-trauma versus trauma patient
self.p_trauma_dist = Bernoulli(
g.prob_trauma, random_seed=self.run_number * g.random_number_set
)
# init sampling for non-stationary poisson process
self.init_nspp()
def init_nspp(self):
# read arrival profile
self.arrivals = pd.read_csv(g.arrival_df) # pylint: disable=attribute-defined-outside-init
self.arrivals["mean_iat"] = 60 / self.arrivals["arrival_rate"]
# maximum arrival rate (smallest time between arrivals)
self.lambda_max = self.arrivals["arrival_rate"].max() # pylint: disable=attribute-defined-outside-init
# thinning exponential
self.arrival_dist = Exponential(
60.0 / self.lambda_max, # pylint: disable=attribute-defined-outside-init
random_seed=self.run_number * g.random_number_set,
)
# thinning uniform rng
self.thinning_rng = Uniform(
low=0.0,
high=1.0, # pylint: disable=attribute-defined-outside-init
random_seed=self.run_number * g.random_number_set,
)
def init_resources(self):
"""
Init the number of resources
and store in the arguments container object
Resource list:
1. Nurses/treatment bays (same thing in this model)
"""
# Shared Resources
self.triage_cubicles = VidigiStore(
self.env, num_resources=g.n_triage, label="triage"
)
self.registration_cubicles = VidigiStore(
self.env, num_resources=g.n_reg, label="registration"
)
# Non-trauma
self.exam_cubicles = VidigiStore(self.env, num_resources=g.n_exam, label="exam")
self.non_trauma_treatment_cubicles = VidigiStore(
self.env, g.n_cubicles_non_trauma_treat, label="non_trauma_treatment"
)
# Trauma
self.trauma_stabilisation_bays = VidigiStore(
self.env, num_resources=g.n_trauma, label="trauma_stabilisation"
)
self.trauma_treatment_cubicles = VidigiStore(
self.env, num_resources=self.n_cubicles_trauma_treat, label="trauma_treatment"
)
# A generator function that represents the DES generator for patient
# arrivals
def generator_patient_arrivals(self):
# We use an infinite loop here to keep doing this indefinitely whilst
# the simulation runs
while True:
t = int(self.env.now // 60) % self.arrivals.shape[0]
lambda_t = self.arrivals["arrival_rate"].iloc[t]
# set to a large number so that at least 1 sample taken!
u = np.inf
interarrival_time = 0.0
# reject samples if u >= lambda_t / lambda_max
while u >= (lambda_t / self.lambda_max):
interarrival_time += self.arrival_dist.sample()
u = self.thinning_rng.sample()
# Freeze this instance of this function in place until the
# inter-arrival time we sampled above has elapsed. Note - time in
# SimPy progresses in "Time Units", which can represent anything
# you like (just make sure you're consistent within the model)
yield self.env.timeout(interarrival_time)
# Increment the patient counter by 1 (this means our first patient
# will have an ID of 1)
self.patient_counter += 1
# Create a new patient - an instance of the Patient Class we
# defined above. Remember, we pass in the ID when creating a
# patient - so here we pass the patient counter to use as the ID.
p = Patient(self.patient_counter)
self.logger.log_arrival(entity_id=p.identifier, pathway="Shared")
# sample if the patient is trauma or non-trauma
trauma = self.p_trauma_dist.sample()
# Tell SimPy to start up the attend_clinic generator function with
# this patient (the generator function that will model the
# patient's journey through the system)
# and store patient in list for later easy access
if trauma:
# create and store a trauma patient to update KPIs.
self.trauma_patients.append(p)
self.env.process(self.attend_trauma_pathway(p))
else:
# create and store a non-trauma patient to update KPIs.
self.non_trauma_patients.append(p)
self.env.process(self.attend_non_trauma_pathway(p))
# A generator function that represents the pathway for a patient going
# through the clinic.
# The patient object is passed in to the generator function so we can
# extract information from / record information to it
def attend_non_trauma_pathway(self, patient):
"""
simulates the non-trauma/minor treatment process for a patient
1. request and wait for sign-in/triage
2. patient registration
3. examination
4a. percentage discharged
4b. remaining percentage treatment then discharge
"""
# record the time of arrival and entered the triage queue
patient.arrival = self.env.now
self.logger.log_queue(
entity_id=patient.identifier,
pathway="Non-Trauma",
event="triage_wait_begins",
)
###################################################
# request sign-in/triage
with self.triage_cubicles.request() as req:
triage_resource = yield req
# record the waiting time for triage
patient.wait_triage = self.env.now - patient.arrival
self.logger.log_resource_use_start(
entity_id=patient.identifier,
pathway="Non-Trauma",
event="triage_begins",
resource_id=triage_resource.id,
unique_resource_id=triage_resource.unique_id,
)
# sample triage duration.
patient.triage_duration = self.triage_dist.sample()
yield self.env.timeout(patient.triage_duration)
self.logger.log_resource_use_end(
entity_id=patient.identifier,
pathway="Non-Trauma",
event="triage_complete",
resource_id=triage_resource.id,
unique_resource_id=triage_resource.unique_id,
)
#########################################################
# record the time that entered the registration queue
start_wait = self.env.now
self.logger.log_queue(
entity_id=patient.identifier,
pathway="Non-Trauma",
event="MINORS_registration_wait_begins",
)
#########################################################
# request registration clerk
with self.registration_cubicles.request() as req:
registration_resource = yield req
# record the waiting time for registration
patient.wait_reg = self.env.now - start_wait
self.logger.log_resource_use_start(
entity_id=patient.identifier,
pathway="Non-Trauma",
event="MINORS_registration_begins",
resource_id=registration_resource.id,
unique_resource_id=registration_resource.unique_id,
)
# sample registration duration.
patient.reg_duration = self.reg_dist.sample()
yield self.env.timeout(patient.reg_duration)
self.logger.log_resource_use_end(
entity_id=patient.identifier,
pathway="Non-Trauma",
event="MINORS_registration_complete",
resource_id=registration_resource.id,
unique_resource_id=registration_resource.unique_id,
)
########################################################
# record the time that entered the evaluation queue
start_wait = self.env.now
self.logger.log_queue(
entity_id=patient.identifier,
pathway="Non-Trauma",
event="MINORS_examination_wait_begins",
)
#########################################################
# request examination resource
with self.exam_cubicles.request() as req:
examination_resource = yield req
# record the waiting time for examination to begin
patient.wait_exam = self.env.now - start_wait
self.logger.log_resource_use_start(
entity_id=patient.identifier,
pathway="Non-Trauma",
event="MINORS_examination_begins",
resource_id=examination_resource.id,
unique_resource_id=examination_resource.unique_id,
)
# sample examination duration.
patient.exam_duration = self.exam_dist.sample()
yield self.env.timeout(patient.exam_duration)
self.logger.log_resource_use_end(
entity_id=patient.identifier,
pathway="Non-Trauma",
event="MINORS_examination_complete",
resource_id=examination_resource.id,
unique_resource_id=examination_resource.unique_id,
)
############################################################################
# sample if patient requires treatment?
patient.require_treat = self.nt_p_treat_dist.sample() # pylint: disable=attribute-defined-outside-init
if patient.require_treat:
# log_custom_event() skips the "unrecognized event_type" warning that
# log_event() would otherwise give for a user-defined event_type like this one.
self.logger.log_custom_event(
entity_id=patient.identifier,
pathway="Non-Trauma",
event="requires_treatment",
event_type="attribute_assigned",
)
# record the time that entered the treatment queue
start_wait = self.env.now
self.logger.log_queue(
entity_id=patient.identifier,
pathway="Non-Trauma",
event="MINORS_treatment_wait_begins",
)
###################################################
# request treatment cubicle
with self.non_trauma_treatment_cubicles.request() as req:
non_trauma_treatment_resource = yield req
# record the waiting time for treatment
patient.wait_treat = self.env.now - start_wait
self.logger.log_resource_use_start(
entity_id=patient.identifier,
pathway="Non-Trauma",
event="MINORS_treatment_begins",
resource_id=non_trauma_treatment_resource.id,
unique_resource_id=non_trauma_treatment_resource.unique_id,
)
# sample treatment duration.
patient.treat_duration = self.nt_treat_dist.sample()
yield self.env.timeout(patient.treat_duration)
self.logger.log_resource_use_end(
entity_id=patient.identifier,
pathway="Non-Trauma",
event="MINORS_treatment_complete",
resource_id=non_trauma_treatment_resource.id,
unique_resource_id=non_trauma_treatment_resource.unique_id,
)
##########################################################################
# Return to what happens to all patients, regardless of whether
# they were sampled as needing treatment
self.logger.log_departure(entity_id=patient.identifier, pathway="Non-Trauma")
# total time in system
patient.total_time = self.env.now - patient.arrival
def attend_trauma_pathway(self, patient):
"""
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
patient.arrival = self.env.now
self.logger.log_queue(
entity_id=patient.identifier, pathway="Trauma", event="triage_wait_begins"
)
###################################################
# request sign-in/triage
with self.triage_cubicles.request() as req:
triage_resource = yield req
# record the waiting time for triage
patient.wait_triage = self.env.now - patient.arrival
self.logger.log_resource_use_start(
entity_id=patient.identifier,
pathway="Trauma",
event="triage_begins",
resource_id=triage_resource.id,
unique_resource_id=triage_resource.unique_id,
)
# sample triage duration.
patient.triage_duration = self.triage_dist.sample()
yield self.env.timeout(patient.triage_duration)
self.logger.log_resource_use_end(
entity_id=patient.identifier,
pathway="Trauma",
event="triage_complete",
resource_id=triage_resource.id,
unique_resource_id=triage_resource.unique_id,
)
###################################################
# record the time that entered the trauma queue
self.logger.log_queue(
entity_id=patient.identifier,
pathway="Trauma",
event="TRAUMA_stabilisation_wait_begins",
)
start_wait = self.env.now
###################################################
# request trauma room
with self.trauma_stabilisation_bays.request() as req:
trauma_resource = yield req
self.logger.log_resource_use_start(
entity_id=patient.identifier,
pathway="Trauma",
event="TRAUMA_stabilisation_begins",
resource_id=trauma_resource.id,
unique_resource_id=trauma_resource.unique_id,
)
# record the waiting time for trauma
patient.wait_trauma = self.env.now - start_wait
# sample stablisation duration.
patient.trauma_duration = self.trauma_dist.sample()
yield self.env.timeout(patient.trauma_duration)
self.logger.log_resource_use_end(
entity_id=patient.identifier,
pathway="Trauma",
event="TRAUMA_stabilisation_complete",
resource_id=trauma_resource.id,
unique_resource_id=trauma_resource.unique_id,
)
#######################################################
# record the time that patient entered the treatment queue
start_wait = self.env.now
self.logger.log_queue(
entity_id=patient.identifier,
pathway="Trauma",
event="TRAUMA_treatment_wait_begins",
)
########################################################
# request treatment cubicle
with self.trauma_treatment_cubicles.request() as req:
trauma_treatment_resource = yield req
# record the waiting time for trauma
patient.wait_treat = self.env.now - start_wait
self.logger.log_resource_use_start(
entity_id=patient.identifier,
pathway="Trauma",
event="TRAUMA_treatment_begins",
resource_id=trauma_treatment_resource.id,
unique_resource_id=trauma_treatment_resource.unique_id,
)
# sample treatment duration.
patient.treat_duration = self.trauma_dist.sample()
yield self.env.timeout(patient.treat_duration)
self.logger.log_resource_use_end(
entity_id=patient.identifier,
pathway="Trauma",
event="TRAUMA_treatment_complete",
resource_id=trauma_treatment_resource.id,
unique_resource_id=trauma_treatment_resource.unique_id,
)
self.logger.log_departure(entity_id=patient.identifier, pathway="Shared")
#########################################################
# total time in system
patient.total_time = self.env.now - patient.arrival
# 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)Show the trial class code
class Trial:
def __init__(self, n_cubicles_trauma_treat=None):
self.n_cubicles_trauma_treat = n_cubicles_trauma_treat
self.all_event_logs = []
self.trial_results_df = pd.DataFrame()
self.run_trial()
# 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(1, g.number_of_runs + 1):
my_model = Model(run, n_cubicles_trauma_treat=self.n_cubicles_trauma_treat)
my_model.run()
self.all_event_logs.append(my_model.logger)
self.trial_results = pd.concat(
[run_results.to_dataframe() for run_results in self.all_event_logs]
)advanced_clinic_simulation = Trial()Constructing a TrialLogger
trial_logs = TrialLogger(
advanced_clinic_simulation.all_event_logs, scenario=g(), label="base scenario"
)
trial_logs<vidigi.logging.TrialLogger at 0x7f9a399b3e30>add_log
A TrialLogger doesn’t have to be built with every run in one go - add_log() appends one more EventLogger to an existing trial. Demonstrated here on a throwaway 5-run copy (rather than the trial_logs used for the rest of this notebook) so the “100 runs” figures quoted further down stay accurate:
demo_trial_logs = TrialLogger(advanced_clinic_simulation.all_event_logs[:5])
extra_run = Model(101)
extra_run.run()
demo_trial_logs.add_log(extra_run.logger)
demo_trial_logs.summary(){'number_of_runs': 6, 'label': None, 'scenario_attached': False}trial_logs.summary(){'number_of_runs': 100, 'label': 'base scenario', 'scenario_attached': True}trial_logs.get_log_by_run(run=2)<vidigi.logging.EventLogger at 0x7f9a3b6826f0>to_dataframe
get_log_by_run above returns one run at a time. to_dataframe() returns every run’s events concatenated into a single DataFrame instead - the same combined frame every plotting/analysis method below builds internally:
trial_logs.to_dataframe()| entity_id | event_type | event | time | pathway | run_number | resource_id | unique_resource_id | |
|---|---|---|---|---|---|---|---|---|
| 0 | 1 | arrival_departure | arrival | 37.593555 | Shared | 1 | NaN | NaN |
| 1 | 1 | queue | triage_wait_begins | 37.593555 | Non-Trauma | 1 | NaN | NaN |
| 2 | 1 | resource_use | triage_begins | 37.593555 | Non-Trauma | 1 | 1.0 | triage_1 |
| 3 | 2 | arrival_departure | arrival | 51.835879 | Shared | 1 | NaN | NaN |
| 4 | 2 | queue | triage_wait_begins | 51.835879 | Non-Trauma | 1 | NaN | NaN |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 146771 | 87 | resource_use_end | MINORS_examination_complete | 598.990148 | Non-Trauma | 100 | 1.0 | exam_1 |
| 146772 | 87 | attribute_assigned | requires_treatment | 598.990148 | Non-Trauma | 100 | NaN | NaN |
| 146773 | 87 | queue | MINORS_treatment_wait_begins | 598.990148 | Non-Trauma | 100 | NaN | NaN |
| 146774 | 87 | resource_use | MINORS_treatment_begins | 598.990148 | Non-Trauma | 100 | 1.0 | non_trauma_treatment_1 |
| 146775 | 86 | resource_use | MINORS_examination_begins | 598.990148 | Non-Trauma | 100 | 1.0 | exam_1 |
146776 rows × 8 columns
Exporting and reloading: to_pickle / read_pickle
to_pickle() pickles the whole TrialLogger - every constituent EventLogger, plus any attached scenario/label - to a file path or a writable binary buffer; read_pickle() loads it back. Handy for caching a long trial’s results, or handing them to another process or notebook without re-running the simulation:
pickle_buffer = io.BytesIO()
trial_logs.to_pickle(pickle_buffer)
pickle_buffer.seek(0)
reloaded_trial_logs = TrialLogger.read_pickle(pickle_buffer)
reloaded_trial_logs.summary(){'number_of_runs': 100, 'label': 'base scenario', 'scenario_attached': True}Feeding the animation
A TrialLogger can build the animation itself - no need to import animate_activity_log or assemble a DataFrame first. run_number= chooses which replication to animate; passing a multi-run TrialLogger without it raises a ValueError that lists the runs available. See the EventLogger notebook for the full animation walkthrough, including event_position_df construction - we reuse the same layout here, since it’s the same model.
event_position_df = create_event_position_df(
[
EventPosition(event="arrival", x=10, y=250, label="Arrival"),
# Triage - minor and trauma
EventPosition(
event="triage_wait_begins", x=160, y=375, label="Waiting for<br>Triage"
),
EventPosition(
event="triage_begins",
x=160,
y=315,
resource="n_triage",
label="Being Triaged",
),
# Minors (non-trauma) pathway
EventPosition(
event="MINORS_registration_wait_begins",
x=300,
y=145,
label="Waiting for<br>Registration",
),
EventPosition(
event="MINORS_registration_begins",
x=300,
y=85,
resource="n_reg",
label="Being<br>Registered",
),
EventPosition(
event="MINORS_examination_wait_begins",
x=465,
y=145,
label="Waiting for<br>Examination",
),
EventPosition(
event="MINORS_examination_begins",
x=465,
y=85,
resource="n_exam",
label="Being<br>Examined",
),
EventPosition(
event="MINORS_treatment_wait_begins",
x=630,
y=145,
label="Waiting for<br>Treatment",
),
EventPosition(
event="MINORS_treatment_begins",
x=630,
y=85,
resource="n_cubicles_non_trauma_treat",
label="Being<br>Treated",
),
# Trauma pathway
EventPosition(
event="TRAUMA_stabilisation_wait_begins",
x=300,
y=560,
label="Waiting for<br>Stabilisation",
),
EventPosition(
event="TRAUMA_stabilisation_begins",
x=300,
y=490,
resource="n_trauma",
label="Being<br>Stabilised",
),
EventPosition(
event="TRAUMA_treatment_wait_begins",
x=630,
y=560,
label="Waiting for<br>Treatment",
),
EventPosition(
event="TRAUMA_treatment_begins",
x=630,
y=490,
resource="n_cubicles_trauma_treat",
label="Being<br>Treated",
),
EventPosition(event="depart", x=670, y=330, label="Exit"),
]
)reshape_for_animations() is available the same way if you want the intermediate frame - the entry point to the step-by-step pipeline - to inspect or tweak before running generate_animation_df()/generate_animation() yourself. It takes the same run_number=:
trial_logs.reshape_for_animations(run_number=2, every_x_time_units=5).head()| snapshot_time | index | entity_id | event_type | event | time | pathway | run_number | resource_id | unique_resource_id | rank | hidden_run_before | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0 |
| 1 | 5 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0 |
| 2 | 10 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0 |
| 3 | 15 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0 |
| 4 | 20 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 0 |
trial_logs.animate_activity_log(
event_position_df,
scenario=g(),
run_number=2,
debug_mode=True,
setup_mode=False,
every_x_time_units=5,
include_play_button=True,
gap_between_entities=11,
gap_between_resources=15,
gap_between_resource_rows=30,
gap_between_queue_rows=30,
plotly_height=600,
plotly_width=1000,
override_x_max=700,
override_y_max=675,
entity_icon_size=10,
resource_icon_size=13,
text_size=15,
wrap_queues_at=10,
step_snapshot_max=20,
limit_duration=g.sim_duration,
time_display_units="dhm",
display_stage_labels=False,
add_background_image="https://raw.githubusercontent.com/Bergam0t/vidigi/refs/heads/main/examples/example_2_branching_multistep/Full%20Model%20Background%20Image%20-%20Horizontal%20Layout.drawio.png",
)Animation function called at 14:07:42
Iteration through time-unit-by-time-unit logs complete 14:07:43
Snapshot df concatenation complete at 14:07:43
Reshaped animation dataframe finished construction at 14:07:43
Placement dataframe started construction at 14:07:43
Placement dataframe finished construction at 14:07:43
Output animation generation complete at 14:07:44
Total Time Elapsed: 1.94 secondsgenerate_dfg
generate_dfg() turns the trial into a directly-follows-graph process map instead of an animation - a diagram overview of the model’s flow. With no run_number=/across_runs= it renders the representative run - the replication whose mean time in system is closest to the trial median. run_number= picks one specific replication, and across_runs=True combines every replication into one map. See the Process Map example for the full walkthrough of all three, plus occupancy_metrics=:
trial_logs.generate_dfg()
TipPrefer the module-level function, or already have a DataFrame?
vidigi.animation.animate_activity_log(event_log=..., ...) still works, and accepts an EventLogger, a TrialLogger, or a plain single-run DataFrame directly:
from vidigi.animation import animate_activity_log
animate_activity_log(
event_log=trial_logs,
event_position_df=event_position_df,
run_number=2,
)The manual route still works too - trial_logs.get_log_by_run(run=2, as_df=True) returns the single-run DataFrame the animation functions expect.
get_event_durations
Every get_event_duration_stat call below summarises this same table - one row per matched pair of events, per entity, per run. Useful on its own for custom analysis get_event_duration_stat doesn’t cover directly:
trial_logs.get_event_durations(
"TRAUMA_treatment_wait_begins", "TRAUMA_treatment_begins"
).head()| entity_id | run_number | pathway | occurrence | first_time | second_time | duration | |
|---|---|---|---|---|---|---|---|
| 0 | 6 | 1 | Trauma | 0 | 357.973070 | 357.973070 | 0.0 |
| 1 | 12 | 1 | Trauma | 0 | 374.484959 | 374.484959 | 0.0 |
| 2 | 23 | 1 | Trauma | 0 | 256.762648 | 256.762648 | 0.0 |
| 3 | 24 | 1 | Trauma | 0 | 445.525327 | 445.525327 | 0.0 |
| 4 | 32 | 1 | Trauma | 0 | 405.569561 | 405.569561 | 0.0 |
trial_logs.get_event_duration_stat(
first_event="TRAUMA_treatment_wait_begins",
second_event="TRAUMA_treatment_begins",
what="mean",
exclude_incomplete=True,
)np.float64(1.45)trial_logs.get_event_duration_stat(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
what="mean",
exclude_incomplete=False,
)np.float64(nan)trial_logs.get_event_duration_stat(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
what="median",
exclude_incomplete=True,
)np.float64(0.0)trial_logs.get_event_duration_stat(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
what="median",
exclude_incomplete=False,
)nantrial_logs.get_event_duration_stat(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
what="count",
exclude_incomplete=True,
)np.int64(1127)trial_logs.get_event_duration_stat(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
what="count",
exclude_incomplete=False,
)1147trial_logs.get_event_duration_stat(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
what="quantile",
exclude_incomplete=True,
q=0.25,
)np.float64(0.0)trial_logs.get_event_duration_stat(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
what="quantile",
exclude_incomplete=True,
q=0.75,
)np.float64(0.0)trial_logs.get_event_duration_stat(
"TRAUMA_treatment_wait_begins", "TRAUMA_treatment_begins", what="served_count"
)np.int64(1127)trial_logs.get_event_duration_stat(
"TRAUMA_treatment_wait_begins", "TRAUMA_treatment_begins", what="summary"
){'mean (of complete)': np.float64(1.45),
'median (of complete)': np.float64(0.0),
'min': np.float64(0.0),
'max': np.float64(94.2),
'unserved_count': np.int64(20),
'served_count': np.int64(1127),
'unserved_rate': np.float64(0.02),
'served_rate': np.float64(0.98),
'unserved_count_mean_per_run': np.float64(0.2),
'served_count_mean_per_run': np.float64(11.27)}get_event_duration_ci
The headline “what is this number, and how sure are we” summary for a trial: the chosen statistic computed separately within each run, then a Student’s t confidence interval taken over those per-replication values - the same computation get_event_duration_stat(..., across="runs") uses for its point estimate, and plot_metric(..., across="runs", error_bars="ci") above draws as an error bar, but returned here as the interval itself rather than a single number or a picture:
trial_logs.get_event_duration_ci(
"TRAUMA_treatment_wait_begins", "TRAUMA_treatment_begins", what="mean"
)ConfidenceInterval(mean=np.float64(1.2082389850423558), half_width=np.float64(0.7875147377275875), lower=np.float64(0.4207242473147683), upper=np.float64(1.9957537227699433), n=100, method='t')get_event_occurrence_rate
A different, per-run question from everything above: in what proportion of replications did an event happen at all, at least once - rather than a per-entity duration or count. "requires_treatment" (the log_custom_event flag logged earlier) turns out to occur in all 100 of 100 runs here - common enough that the point estimate is a clean 100%, though the Wilson interval’s lower bound (96.3%) is still shy of certainty from only 100 observations. A genuinely rare, binary condition (a capacity breach, a specific alarm) would show the opposite shape: a low point estimate with a proportionally wide interval, since Wilson needs to stay inside [0, 1] and behave sensibly near the edges - unlike a Student’s t interval, which could stray outside [0, 1] for a rate this close to 100%:
trial_logs.get_event_occurrence_rate("requires_treatment")ProportionEstimate(proportion=1.0, lower=np.float64(0.9630065017930143), upper=1.0, n_runs=100, n_occurred=100, ci_level=0.95, method='wilson')event_pairs = [
{
"first_event": "triage_wait_begins",
"second_event": "triage_begins",
"label": "Triage Wait Length",
},
{
"first_event": "MINORS_registration_wait_begins",
"second_event": "MINORS_registration_begins",
"label": "Minors Registration Wait Length",
},
{
"first_event": "MINORS_examination_wait_begins",
"second_event": "MINORS_examination_begins",
"label": "Minors Examination Wait Length",
},
{
"first_event": "MINORS_treatment_wait_begins",
"second_event": "MINORS_treatment_begins",
"label": "Minors Treatment Wait Length",
},
{
"first_event": "TRAUMA_stabilisation_wait_begins",
"second_event": "TRAUMA_stabilisation_begins",
"label": "Trauma Stabilisation Wait Length",
},
{
"first_event": "TRAUMA_treatment_wait_begins",
"second_event": "TRAUMA_treatment_begins",
"label": "Trauma Treatment Wait Length",
},
]
fig = trial_logs.plot_metric(event_pairs, kind="bar", what="mean")
fig.update_layout(title="Mean step durations", width=800)fig = trial_logs.plot_metric(event_pairs, kind="bar", what="median")
fig.update_layout(title="Median step durations", width=800)fig = trial_logs.plot_metric(event_pairs, kind="bar", what="max")
fig.update_layout(title="Max step durations", width=800)trial_logs.plot_queue_size(
["triage_wait_begins"],
limit_duration=g.sim_duration,
every_x_time_units=30,
width=1000,
title="Triage Queue Length",
)trial_logs.plot_queue_size(
[
"triage_wait_begins",
"MINORS_registration_wait_begins",
"MINORS_examination_wait_begins",
"MINORS_examination_begins",
"MINORS_treatment_wait_begins",
"TRAUMA_stabilisation_wait_begins",
"TRAUMA_treatment_wait_begins",
],
limit_duration=g.sim_duration,
every_x_time_units=30,
show_all_runs=False,
height=1200,
width=1200,
)/home/runner/work/vidigi/vidigi/src/vidigi/plots.py:205: UserWarning: ['MINORS_treatment_wait_begins'] did not occur in any run, so will be reported as a queue of zero throughout. Check these against the values in your event log's 'event' column.
event_counts = queue_size_over_time(highlight_bands (new in 2.0.0) shades a threshold zone behind the chart - here, a red zone flagging a queue of 10 or more. shared_y_axis=True (the default) means every facet uses the same y-range, so one band placed at the shared scale reads consistently across all seven steps, and is drawn on every facet rather than only the first:
trial_logs.plot_queue_size(
[
"triage_wait_begins",
"MINORS_registration_wait_begins",
"MINORS_examination_wait_begins",
"MINORS_examination_begins",
"MINORS_treatment_wait_begins",
"TRAUMA_stabilisation_wait_begins",
"TRAUMA_treatment_wait_begins",
],
limit_duration=g.sim_duration,
every_x_time_units=30,
show_all_runs=True,
height=1200,
width=1200,
highlight_bands=[{"lower": 10, "colour": "red", "label": "long queue (>=10)"}],
)/home/runner/work/vidigi/vidigi/src/vidigi/plots.py:205: UserWarning: ['MINORS_treatment_wait_begins'] did not occur in any run, so will be reported as a queue of zero throughout. Check these against the values in your event log's 'event' column.
event_counts = queue_size_over_time(
TipDeciding how much of that queue is warm-up
The queue lengths above include the startup transient - the artificially short queue every run starts with, because nobody was there before the first arrival. TrialLogger.plot_warm_up_diagnostic() gives a visual way to decide how much of the start of a run to discard before trusting any statistic computed from it. See feat_warm_up.ipynb for the full walkthrough.
Distributions of durations
plot_metric and get_event_duration_stat above only ever give us a single number - the mean, median, or another summary statistic of a duration. For a queueing model that can hide exactly the detail that matters: two systems with the same mean wait can look very different once you see the spread, or whether a handful of patients are waiting far longer than everyone else.
TrialLogger.plot_duration_distribution plots the full distribution of a duration instead, in one of four styles (kind="hist"|"box"|"violin"|"ecdf"). It’s built directly on the same get_event_durations used above, so it never computes anything get_event_duration_stat couldn’t already tell you - it just shows the whole shape rather than reducing it to one number.
trial_logs.plot_duration_distribution(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
kind="hist",
title="Trauma treatment wait - distribution of durations",
)A box or violin plot shows the same distribution more compactly - handy when comparing several durations side by side, which we’ll do below.
trial_logs.plot_duration_distribution(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
kind="box",
title="Trauma treatment wait - box plot",
)trial_logs.plot_duration_distribution(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
kind="violin",
title="Trauma treatment wait - violin plot",
)An empirical cumulative distribution function (ECDF) is useful when you care about a specific threshold - for example, “what proportion of patients waited more than 30 minutes for trauma treatment?” can be read straight off the chart.
trial_logs.plot_duration_distribution(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
kind="ecdf",
title="Trauma treatment wait - ECDF",
)Passing split_by="run" draws one trace per replication instead of pooling every run’s durations together - useful for seeing how much a duration’s spread itself varies from run to run, not just its mean. With 100 replications in this trial that would be a very busy chart, so here we build a small TrialLogger from just the first twenty runs to demonstrate it.
first_five_runs = TrialLogger(advanced_clinic_simulation.all_event_logs[:20])
first_five_runs.plot_duration_distribution(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
kind="box",
split_by="run",
title="Trauma treatment wait by run (first 20 runs)",
)That’s a reasonable way to look at a handful of runs, but it doesn’t scale to all 100 - 100 overlapping box plots would be unreadable. Two more kind options are built for exactly this case:
kind="ridgeline"stacks one density curve per group with a slight vertical overlap (a “joy plot”). It reads well up to a few dozen groups, but its offset-stacking is linear in the number of groups - stretched over all 100 runs the curves and y-axis labels overlap into an unreadable smear, so we draw it against the same 20-run subset as the box plot above rather than the full trial.kind="heatmap"draws duration along the x-axis and one row per group, coloured by density - since it costs no vertical space per row at all, this is the one that actually scales to the full 100 runs, and is the better choice whenever a ridgeline would get too tall to read.
fig = first_five_runs.plot_duration_distribution(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
kind="ridgeline",
split_by="run",
title="Trauma treatment wait by run - ridgeline (first 20 runs)",
)
fig.update_layout(height=800)fig = trial_logs.plot_duration_distribution(
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
kind="heatmap",
split_by="run",
title="Trauma treatment wait by run - heatmap (all 100 runs)",
)
fig.update_layout(height=600)Uncertainty on a bar chart
Every plot_metric call above pools every entity’s duration into one statistic per bar - across="entities", the default, matching every release before 2.0.0. For a stochastic model that hides how much a step’s duration actually varies between replications, which is what determines how much you should trust the bar height.
across="runs" computes the statistic separately within each run first, then draws the mean of those per-run values - the number an error_bars interval is actually about. Entities within a run are correlated (a bad morning makes many waits long together), so an interval would be badly overconfident if it pretended every entity was an independent observation; replications are the independent unit here, not entities. error_bars="ci" needs the optional scipy dependency (pip install vidigi[stats]) - "sd", "se" and the asymmetric "range"/"iqr" do not.
fig = trial_logs.plot_metric(
event_pairs,
kind="bar",
what="mean",
across="runs",
error_bars="ci",
show_runs=True,
)
fig.update_layout(title="Mean step durations, with 95% CI across runs", width=900)The scattered points are each run’s own mean, and the error bar is a 95% confidence interval computed over them (ci_level= changes the level). Compare this with the pooled-entity bar chart drawn near the top of this notebook - here, a step with high run-to-run variability gets a visibly wide interval even where its pooled mean alone would look no different from a more consistent step.
Spotting a fluke replication: get_outlier_runs / plot_outlier_runs
More replications means more chances one of them is a fluke - a run whose per-replication value happens to land well away from the rest, by chance rather than because anything about the model changed. get_outlier_runs takes the same per-replication values get_event_duration_ci/plot_replication_analysis are built on, and applies Tukey’s classic 1.5x IQR fence - the same convention error_bars="iqr" on plot_metric/plot_resource_utilisation already draws as an error bar, just applied here as a threshold rather than a picture. It returns the whole per-run table, not just a list of flagged run numbers, so the fence values behind each flag are there to inspect too:
trial_logs.get_outlier_runs("MINORS_examination_wait_begins", "MINORS_examination_begins")| run_number | value | lower_fence | upper_fence | is_outlier | |
|---|---|---|---|---|---|
| 0 | 1 | 24.744272 | 10.52118 | 42.400734 | False |
| 1 | 2 | 15.222273 | 10.52118 | 42.400734 | False |
| 2 | 3 | 16.578365 | 10.52118 | 42.400734 | False |
| 3 | 4 | 25.097810 | 10.52118 | 42.400734 | False |
| 4 | 5 | 23.123154 | 10.52118 | 42.400734 | False |
| ... | ... | ... | ... | ... | ... |
| 95 | 96 | 29.942848 | 10.52118 | 42.400734 | False |
| 96 | 97 | 28.388232 | 10.52118 | 42.400734 | False |
| 97 | 98 | 27.263757 | 10.52118 | 42.400734 | False |
| 98 | 99 | 33.381774 | 10.52118 | 42.400734 | False |
| 99 | 100 | 27.460647 | 10.52118 | 42.400734 | False |
100 rows × 5 columns
plot_outlier_runs draws the same fence as a picture instead of a table: a horizontal beeswarm of every run’s value, red-shaded bands marking the zones beyond the lower/upper fence, and each point coloured (and shape-coded, for a colourblind-safe read) by whether it was flagged. The title states the same verdict the table above gives, as a sentence:
fig = trial_logs.plot_outlier_runs(
"MINORS_examination_wait_begins", "MINORS_examination_begins"
)
fig.update_layout(width=900, height=400)
TipTwo more diagnostics that live in their own notebooks
TrialLogger.plot_replication_analysis()/get_replication_precision() ask a related question to the outlier check above - not “is one run a fluke” but “have I run enough replications at all” - watching a confidence interval’s width shrink as runs accumulate. See feat_replication_analysis.ipynb for the full walkthrough.
TrialLogger.plot_metric_vs_arrival_time()/get_entity_metric_by_arrival() ask yet another: within one run’s steady operation, does a metric drift depending on when the entity that produced it arrived. See feat_metric_vs_arrival_time.ipynb.
Resource utilisation
resource_use/resource_use_end events are already in the log - every with self.*.request() block above logs one pair - but computing utilisation from them needs one more thing the log does not carry: how many units of each resource exist. get_resource_utilisation resolves that from a scenario object plus a resource_map naming which attribute holds each step’s capacity, reusing the same idea as the animation’s resource= argument on EventPosition. (resource_capacities={step: count} works too, if you don’t have a scenario object handy; see vidigi.analysis._resolve_resource_capacities for the full set of options.)
by="step" (the default) gives one row per run per step:
resource_map = {
"triage_begins": "n_triage",
"MINORS_registration_begins": "n_reg",
"MINORS_examination_begins": "n_exam",
"MINORS_treatment_begins": "n_cubicles_non_trauma_treat",
"TRAUMA_stabilisation_begins": "n_trauma",
"TRAUMA_treatment_begins": "n_cubicles_trauma_treat",
}
utilisation_by_step = trial_logs.get_resource_utilisation(
by="step", scenario=g(), resource_map=resource_map, limit_duration=g.sim_duration
)
# One row per run per step - summarise across runs to see the shape of the result
utilisation_by_step.groupby("event")[["mean_in_use", "utilisation"]].mean()/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 1398 resource_use row(s) were still open at the end of the window - censored at 600 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(| mean_in_use | utilisation | |
|---|---|---|
| event | ||
| MINORS_examination_begins | 2.139331 | 0.713110 |
| MINORS_registration_begins | 1.336405 | 0.668202 |
| MINORS_treatment_begins | 1.035854 | 0.258964 |
| TRAUMA_stabilisation_begins | 1.751918 | 0.437979 |
| TRAUMA_treatment_begins | 1.220636 | 0.244127 |
| triage_begins | 1.286545 | 0.643273 |
mean_in_use needs no capacity at all (it’s just busy time divided by the window length) and is always populated; utilisation additionally divides by the resolved capacity. by="resource" instead gives one row per run per physical unit, where capacity is always exactly 1 - no capacity route is needed at all:
trial_logs.get_resource_utilisation(by="resource", limit_duration=g.sim_duration).head()/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 1398 resource_use row(s) were still open at the end of the window - censored at 600 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(| run_number | resource_id | busy_time | mean_in_use | capacity | utilisation | |
|---|---|---|---|---|---|---|
| 0 | 1 | exam_1 | 458.954452 | 0.764924 | 1.0 | 0.764924 |
| 1 | 1 | exam_2 | 444.995694 | 0.741659 | 1.0 | 0.741659 |
| 2 | 1 | exam_3 | 445.538956 | 0.742565 | 1.0 | 0.742565 |
| 3 | 1 | non_trauma_treatment_1 | 157.403616 | 0.262339 | 1.0 | 0.262339 |
| 4 | 1 | non_trauma_treatment_2 | 152.961592 | 0.254936 | 1.0 | 0.254936 |
Notice several of those utilisation values are above 1 - not a genuine over-capacity reading, but the resource_id collision this model would have if it only logged that column. by="resource" assumes resource_id is unique across the whole log, but this model’s VidigiStore instances each number their own units starting from 1 - triage cubicle 1 and registration clerk 1 share the same resource_id, so their busy time gets silently summed as if they were one physical thing.
The model code above avoids this: every VidigiStore is constructed with a label= (e.g. label="triage"), and every log_resource_use_start/log_resource_use_end call also logs unique_resource_id=<resource>.unique_id alongside the animation-safe resource_id. Pointing by="resource" at that column instead - via resource_col_name="unique_resource_id" - resolves the collision:
trial_logs.get_resource_utilisation(
by="resource", resource_col_name="unique_resource_id", limit_duration=g.sim_duration
).head()/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 1398 resource_use row(s) were still open at the end of the window - censored at 600 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(| run_number | resource_id | busy_time | mean_in_use | capacity | utilisation | |
|---|---|---|---|---|---|---|
| 0 | 1 | exam_1 | 458.954452 | 0.764924 | 1.0 | 0.764924 |
| 1 | 1 | exam_2 | 444.995694 | 0.741659 | 1.0 | 0.741659 |
| 2 | 1 | exam_3 | 445.538956 | 0.742565 | 1.0 | 0.742565 |
| 3 | 1 | non_trauma_treatment_1 | 157.403616 | 0.262339 | 1.0 | 0.262339 |
| 4 | 1 | non_trauma_treatment_2 | 152.961592 | 0.254936 | 1.0 | 0.254936 |
by="run" pools every step into a single blended figure per run. This model has six genuinely different resource types (triage nurses, registration clerks, exam rooms, trauma bays, and two kinds of treatment cubicle), so pooling them warns - a single “how utilised was everything” number is rarely what a capacity-planning question is actually asking when the resources aren’t interchangeable:
trial_logs.get_resource_utilisation(
by="run", scenario=g(), resource_map=resource_map, limit_duration=g.sim_duration
).head()/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 1398 resource_use row(s) were still open at the end of the window - censored at 600 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(
/home/runner/work/vidigi/vidigi/src/vidigi/logging.py:1805: UserWarning: by='run' is pooling 6 distinct steps (['MINORS_examination_begins', 'MINORS_registration_begins', 'MINORS_treatment_begins', 'TRAUMA_stabilisation_begins', 'TRAUMA_treatment_begins', 'triage_begins']) into one blended busy-time and utilisation figure. This treats every unit of every step as interchangeable, which is rarely the number a capacity-planning question is asking for a system with more than one resource type - consider by='step' instead.
return resource_utilisation(| run_number | busy_time | mean_in_use | capacity | utilisation | |
|---|---|---|---|---|---|
| 0 | 1 | 5193.255354 | 8.655426 | 20 | 0.432771 |
| 1 | 2 | 5123.119608 | 8.538533 | 20 | 0.426927 |
| 2 | 3 | 4692.159652 | 7.820266 | 20 | 0.391013 |
| 3 | 4 | 5413.127194 | 9.021879 | 20 | 0.451094 |
| 4 | 5 | 5558.386619 | 9.263978 | 20 | 0.463199 |
Plotting resource utilisation
plot_resource_utilisation is a thin wrapper over get_resource_utilisation above: the bar height is the mean across runs, and the error bar is a confidence interval over the same per-run values (error_bars="ci" and show_runs=True are the defaults here, unlike plot_metric, since this function is new rather than extracted from older, bar-only behaviour). The dashed line at 1.0 is diagnostic - utilisation can never legitimately exceed it, so a bar crossing it means the resolved capacity, the logged intervals, or both need a second look.
trial_logs.plot_resource_utilisation(
by="step", scenario=g(), resource_map=resource_map, limit_duration=g.sim_duration
)/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 1398 resource_use row(s) were still open at the end of the window - censored at 600 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(by="resource" draws the same chart broken down to one bar per physical unit instead of per step. Using the default resource_id column, several bars cross the dashed line at 1.0 here - the same resource_id collision as above, not a genuine over-capacity reading:
trial_logs.plot_resource_utilisation(by="resource", limit_duration=g.sim_duration)/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 1398 resource_use row(s) were still open at the end of the window - censored at 600 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(Passing resource_col_name="unique_resource_id" through resolves it here too:
trial_logs.plot_resource_utilisation(
by="resource", resource_col_name="unique_resource_id", limit_duration=g.sim_duration
)/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 1398 resource_use row(s) were still open at the end of the window - censored at 600 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(plot_resource_utilisation_over_time is the resource equivalent of plot_queue_size above - how many units of each step were busy at each snapshot, across every run, faceted one panel per step. Traces use a step-function (hv) line shape rather than a straight one, since linear interpolation between snapshots would draw fractional resource counts that never existed.
fig = trial_logs.plot_resource_utilisation_over_time(
every_x_time_units=30, limit_duration=g.sim_duration, show_all_runs=False
)
fig.update_layout(height=1400)/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:2165: UserWarning: 1398 resource_use row(s) were still open at the end of the window - censored at 600 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(Passing as_proportion=True divides each step’s count by its resolved capacity, so every panel shares the same 0-1 scale regardless of how many units that step actually has - handy for comparing a 2-unit resource against a 5-unit one at a glance. Unlike plot_resource_utilisation, there is no fallback here if a plotted step’s capacity can’t be resolved: a partially-NaN proportion trace would be more misleading than an error naming the problem, so it raises instead. highlight_bands works here too, and since every step already shares the same 0-1 proportion scale, one band reads consistently across all six facets:
fig = trial_logs.plot_resource_utilisation_over_time(
every_x_time_units=30,
limit_duration=g.sim_duration,
show_all_runs=False,
as_proportion=True,
scenario=g(),
resource_map=resource_map,
highlight_bands=[{"lower": 0.9, "colour": "red", "label": "near capacity (>=90%)"}],
)
fig.update_layout(height=1400)/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:2165: UserWarning: 1398 resource_use row(s) were still open at the end of the window - censored at 600 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(
/home/runner/work/vidigi/vidigi/src/vidigi/plots.py:1811: UserWarning: 1398 resource_use row(s) were still open at the end of the window - censored at 600 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(Comparing two scenarios
Everything above describes one scenario. The question a capacity-planning exercise usually asks next is comparative: if the clinic ran with fewer trauma treatment cubicles, would patients notice? compare_event_duration_stat - and its resource-utilisation twin, compare_resource_utilisation - take a second TrialLogger and answer exactly that.
Model/Trial above take an optional n_cubicles_trauma_treat= override (a verified no-op when omitted - every run so far used the default g.n_cubicles_trauma_treat, 5). A second trial with 3 cubicles instead of 5 gives a deliberately worse scenario to compare against. Each run reuses the same random seeds as its counterpart in trial_logs, since the model’s distributions are seeded from run_number * g.random_number_set, independent of n_cubicles_trauma_treat - so this is a matched (common-random-numbers) comparison. The reduced scenario gets its own small subclass of g, rather than mutating g.n_cubicles_trauma_treat directly, so trial_logs.scenario.n_cubicles_trauma_treat (attached to the first trial) keeps reading 5:
class GFewerTraumaCubicles(g):
n_cubicles_trauma_treat = 3
reduced_trial = Trial(n_cubicles_trauma_treat=GFewerTraumaCubicles.n_cubicles_trauma_treat)
reduced_trial_logs = TrialLogger(
reduced_trial.all_event_logs,
scenario=GFewerTraumaCubicles(),
label="3 trauma treatment cubicles",
)
reduced_trial_logs.summary(){'number_of_runs': 100,
'label': '3 trauma treatment cubicles',
'scenario_attached': True}compare_event_duration_stat computes a confidence interval on each side independently (never pooled - the two trials are different scenarios, not before/after pairs of the same run), plus a Welch’s t-test p-value alongside it. ci_overlap=False is a safe “these two scenarios differ” signal; ci_overlap=True only means “not conclusively different by this simple check”, not proof they’re the same - p_value is the more rigorous figure to read alongside it. label_a/label_b default to each trial’s own .label, which is why neither is passed below:
comparison = trial_logs.compare_event_duration_stat(
reduced_trial_logs, "TRAUMA_treatment_wait_begins", "TRAUMA_treatment_begins"
)
print(f"{comparison.label_a} mean wait: {comparison.mean_a:.1f}")
print(f"{comparison.label_b} mean wait: {comparison.mean_b:.1f}")
print(f"delta: {comparison.delta:+.1f} ({comparison.delta_pct:+.0f}%)")
print(f"CIs overlap: {comparison.ci_overlap} Welch's t-test p-value: {comparison.p_value:.4g}")base scenario mean wait: 1.2
3 trauma treatment cubicles mean wait: 10.0
delta: +8.8 (+727%)
CIs overlap: False Welch's t-test p-value: 1.675e-07plot_event_duration_comparison draws the same numbers as a two-bar chart, with a CI error bar on each and a title stating the overlap verdict and p-value - worded to avoid overclaiming (“CIs overlap - not conclusively different”, never “no difference”). highlight_bands works here too:
fig = trial_logs.plot_event_duration_comparison(
reduced_trial_logs,
"TRAUMA_treatment_wait_begins",
"TRAUMA_treatment_begins",
highlight_bands=[{"upper": 30, "colour": "green", "label": "acceptable (<30 min)"}],
)
fig.update_layout(width=700, height=500)compare_resource_utilisation is the same idea for utilisation, always pooled across every step/resource into one blended per-run figure (by="run"); call get_resource_utilisation(by=...) on each trial directly and pass the result into vidigi.analysis.compare_replication_values to compare one specific step or resource instead. Each trial resolves capacity from its own attached scenario (g() for trial_logs, GFewerTraumaCubicles() for reduced_trial_logs, both attached at construction above), so only resource_map needs passing here:
util_comparison = trial_logs.compare_resource_utilisation(
reduced_trial_logs, resource_map=resource_map
)
print(f"{util_comparison.label_a} utilisation: {util_comparison.mean_a:.0%}")
print(f"{util_comparison.label_b} utilisation: {util_comparison.mean_b:.0%}")
print(
f"CIs overlap: {util_comparison.ci_overlap} "
f"Welch's t-test p-value: {util_comparison.p_value:.4g}"
)/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 1398 resource_use row(s) were still open at the end of the window - censored at 599.9926563616847 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(
/home/runner/work/vidigi/vidigi/src/vidigi/logging.py:1805: UserWarning: by='run' is pooling 6 distinct steps (['MINORS_examination_begins', 'MINORS_registration_begins', 'MINORS_treatment_begins', 'TRAUMA_stabilisation_begins', 'TRAUMA_treatment_begins', 'triage_begins']) into one blended busy-time and utilisation figure. This treats every unit of every step as interchangeable, which is rarely the number a capacity-planning question is asking for a system with more than one resource type - consider by='step' instead.
return resource_utilisation(base scenario utilisation: 44%
3 trauma treatment cubicles utilisation: 48%
CIs overlap: False Welch's t-test p-value: 1.56e-10/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 1357 resource_use row(s) were still open at the end of the window - censored at 599.9926563616847 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(
/home/runner/work/vidigi/vidigi/src/vidigi/logging.py:1805: UserWarning: by='run' is pooling 6 distinct steps (['MINORS_examination_begins', 'MINORS_registration_begins', 'MINORS_treatment_begins', 'TRAUMA_stabilisation_begins', 'TRAUMA_treatment_begins', 'triage_begins']) into one blended busy-time and utilisation figure. This treats every unit of every step as interchangeable, which is rarely the number a capacity-planning question is asking for a system with more than one resource type - consider by='step' instead.
return resource_utilisation(plot_resource_utilisation_comparison is plot_event_duration_comparison’s twin for utilisation - the same two-bar-plus-CI chart, built on compare_resource_utilisation instead of compare_event_duration_stat, so the title states the same kind of overlap verdict and p-value. highlight_bands works here too, alongside the dashed 100% reference line already drawn on every utilisation chart:
fig = trial_logs.plot_resource_utilisation_comparison(
reduced_trial_logs,
resource_map=resource_map,
highlight_bands=[{"lower": 0.9, "colour": "red", "label": "near capacity (>=90%)"}],
)
fig.update_layout(width=700, height=500)/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 1398 resource_use row(s) were still open at the end of the window - censored at 599.9926563616847 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(
/home/runner/work/vidigi/vidigi/src/vidigi/plots.py:2908: UserWarning: by='run' is pooling 6 distinct steps (['MINORS_examination_begins', 'MINORS_registration_begins', 'MINORS_treatment_begins', 'TRAUMA_stabilisation_begins', 'TRAUMA_treatment_begins', 'triage_begins']) into one blended busy-time and utilisation figure. This treats every unit of every step as interchangeable, which is rarely the number a capacity-planning question is asking for a system with more than one resource type - consider by='step' instead.
values_a = resource_utilisation(
/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 1357 resource_use row(s) were still open at the end of the window - censored at 599.9926563616847 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(
/home/runner/work/vidigi/vidigi/src/vidigi/plots.py:2911: UserWarning: by='run' is pooling 6 distinct steps (['MINORS_examination_begins', 'MINORS_registration_begins', 'MINORS_treatment_begins', 'TRAUMA_stabilisation_begins', 'TRAUMA_treatment_begins', 'triage_begins']) into one blended busy-time and utilisation figure. This treats every unit of every step as interchangeable, which is rarely the number a capacity-planning question is asking for a system with more than one resource type - consider by='step' instead.
values_b = resource_utilisation(That’s the same question compare_event_duration_stat answered above, now for resource utilisation instead of wait time - turning “would patients notice?” into a number, a confidence interval, and a plain-language verdict, rather than eyeballing two separate plot_metric/plot_resource_utilisation charts and guessing.