import plotly.io as pio
from feat_animation_warm_up_model_classes import Trial, g
from vidigi.animation import animate_activity_log
from vidigi.prep import reshape_for_animations
from vidigi.utils import EventPosition, create_event_position_df
pio.renderers.default = "notebook"Feature Example: Discarding a Warm-up Period from an Animation
event_durations, queue_size_over_time and the rest of vidigi.analysis all take a warm_up= parameter for excluding early, unrepresentative results - see feat_warm_up.ipynb for choosing how much. This notebook is about a related but distinct problem: applying that same warm-up to an animation, where the obvious approach - filtering the event log by time - does not just exclude early data, it actively corrupts the result.
Model setup
The same single-step clinic model used by example_1_simplest_case - one queue, one resource, patients arrive, wait for a treatment cubicle, are treated, and leave.
NoteView Imported Code, which has had logging steps added at the appropriate points in the ‘model’ class
import random
import numpy as np
import pandas as pd
import simpy
from sim_tools.distributions import Exponential, Lognormal
from vidigi.resources import VidigiStore
# 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_cubicles: int
The number of treatment cubicles
trauma_treat_mean: float
Mean of the trauma cubicle treatment distribution (Lognormal)
trauma_treat_var: float
Variance of the trauma cubicle treatment distribution (Lognormal)
arrival_rate: float
Set the mean of the exponential distribution that is used to sample the
inter-arrival time of patients
sim_duration: int
The number of time units the simulation will run for
number_of_runs: int
The number of times the simulation will be run with different random number streams
"""
random_number_set = 42
n_cubicles = 4
trauma_treat_mean = 40
trauma_treat_var = 5
arrival_rate = 5
sim_duration = 600
number_of_runs = 100
# Class representing patients coming in to the clinic.
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
self.arrival = -np.inf
self.wait_treat = -np.inf
self.total_time = -np.inf
self.treat_duration = -np.inf
# 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):
# Create a SimPy environment in which everything will live
self.env = simpy.Environment()
self.event_log = []
# Create a patient counter (which we'll use as a patient ID)
self.patient_counter = 0
self.patients = []
# Create our resources
self.init_resources()
# Store the passed in run number
self.run_number = run_number
# Create a new Pandas DataFrame that will store some results against
# the patient ID (which we'll use as the index).
self.results_df = pd.DataFrame()
self.results_df["Patient ID"] = [1]
self.results_df["Queue Time Cubicle"] = [0.0]
self.results_df["Time with Nurse"] = [0.0]
self.results_df.set_index("Patient ID", inplace=True)
# Create an attribute to store the mean queuing times across this run of
# the model
self.mean_q_time_cubicle = 0
self.patient_inter_arrival_dist = Exponential(
mean=g.arrival_rate, random_seed=self.run_number * g.random_number_set
)
self.treat_dist = Lognormal(
mean=g.trauma_treat_mean,
stdev=g.trauma_treat_var,
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)
"""
self.treatment_cubicles = VidigiStore(
self.env, num_resources=g.n_cubicles, label="treatment_cubicle"
)
# 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:
# 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)
# Store patient in list for later easy access
self.patients.append(p)
# 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)
self.env.process(self.attend_clinic(p))
# Randomly sample the time to the next patient arriving. Here, we
# sample from an exponential distribution (common for inter-arrival
# times), and pass in a lambda value of 1 / mean. The mean
# inter-arrival time is stored in the g class.
sampled_inter = self.patient_inter_arrival_dist.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(sampled_inter)
# 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_clinic(self, patient):
self.arrival = self.env.now
self.event_log.append(
{
"patient": patient.identifier,
"pathway": "Simplest",
"event_type": "arrival_departure",
"event": "arrival",
"time": self.env.now,
}
)
# request examination resource
start_wait = self.env.now
self.event_log.append(
{
"patient": patient.identifier,
"pathway": "Simplest",
"event": "treatment_wait_begins",
"event_type": "queue",
"time": self.env.now,
}
)
# Seize a treatment resource when available
treatment_resource = yield self.treatment_cubicles.get_direct()
# record the waiting time for registration
self.wait_treat = self.env.now - start_wait
self.event_log.append(
{
"patient": patient.identifier,
"pathway": "Simplest",
"event": "treatment_begins",
"event_type": "resource_use",
"time": self.env.now,
"resource_id": treatment_resource.id_attribute,
}
)
# sample treatment duration
self.treat_duration = self.treat_dist.sample()
yield self.env.timeout(self.treat_duration)
self.event_log.append(
{
"patient": patient.identifier,
"pathway": "Simplest",
"event": "treatment_complete",
"event_type": "resource_use_end",
"time": self.env.now,
"resource_id": treatment_resource.id_attribute,
}
)
# Resource is no longer in use, so put it back in
self.treatment_cubicles.put(treatment_resource)
# total time in system
self.total_time = self.env.now - self.arrival
self.event_log.append(
{
"patient": patient.identifier,
"pathway": "Simplest",
"event": "depart",
"event_type": "arrival_departure",
"time": self.env.now,
}
)
# This method calculates results over a single run. Here we just calculate
# a mean, but in real world models you'd probably want to calculate more.
def calculate_run_results(self):
# Take the mean of the queuing times across patients in this run of the
# model.
self.mean_q_time_cubicle = self.results_df["Queue Time Cubicle"].mean()
# The run method starts up the DES entity generators, runs the simulation,
# and in turns calls anything we need to generate results for the run
def run(self):
# Start up our DES entity generators that create new patients. We've
# only got one in this model, but we'd need to do this for each one if
# we had multiple generators.
self.env.process(self.generator_patient_arrivals())
# Run the model for the duration specified in g class
self.env.run(until=g.sim_duration)
# Now the simulation run has finished, call the method that calculates
# run results
self.calculate_run_results()
self.event_log = pd.DataFrame(self.event_log)
self.event_log["run"] = self.run_number
return {"results": self.results_df, "event_log": self.event_log}
# Class representing a Trial for our simulation - a batch of simulation runs.
class Trial:
# The constructor sets up a pandas dataframe that will store the key
# results from each run against run number, with run number as the index.
def __init__(self):
self.df_trial_results = pd.DataFrame()
self.df_trial_results["Run Number"] = [0]
self.df_trial_results["Arrivals"] = [0]
self.df_trial_results["Mean Queue Time Cubicle"] = [0.0]
self.df_trial_results.set_index("Run Number", inplace=True)
self.all_event_logs = []
# Method to run a trial
def run_trial(self):
print(f"{g.n_cubicles} nurses")
print() ## Print a blank line
# Run the simulation for the number of runs specified in g class.
# For each run, we create a new instance of the Model class and call its
# run method, which sets everything else in motion. Once the run has
# completed, we grab out the stored run results (just mean queuing time
# here) and store it against the run number in the trial results
# dataframe.
for run in range(g.number_of_runs):
random.seed(run)
my_model = Model(run)
model_outputs = my_model.run()
patient_level_results = model_outputs["results"]
event_log = model_outputs["event_log"]
self.df_trial_results.loc[run] = [
len(patient_level_results),
my_model.mean_q_time_cubicle,
]
# print(event_log)
self.all_event_logs.append(event_log)
self.all_event_logs = pd.concat(self.all_event_logs)g.number_of_runs = 1
my_trial = Trial()
my_trial.run_trial()
event_log = my_trial.all_event_logs[my_trial.all_event_logs["run"] == 0].reset_index(
drop=True
)
event_log.head()4 nurses
| patient | pathway | event_type | event | time | resource_id | run | |
|---|---|---|---|---|---|---|---|
| 0 | 1 | Simplest | arrival_departure | arrival | 0.00000 | NaN | 0 |
| 1 | 1 | Simplest | queue | treatment_wait_begins | 0.00000 | NaN | 0 |
| 2 | 1 | Simplest | resource_use | treatment_begins | 0.00000 | 1.0 | 0 |
| 3 | 2 | Simplest | arrival_departure | arrival | 3.39966 | NaN | 0 |
| 4 | 2 | Simplest | queue | treatment_wait_begins | 3.39966 | NaN | 0 |
4 cubicles, a fairly busy arrival rate, and a 600 time-unit run - a queue is already well established by the time we cut it off at warm_up = 100.
warm_up = 100
still_present_at_warm_up = (
event_log[event_log["event"] == "arrival"][["patient", "time"]]
.rename(columns={"time": "arrival_time"})
.merge(
event_log[event_log["event"] == "depart"][["patient", "time"]].rename(
columns={"time": "depart_time"}
),
on="patient",
how="left",
)
)
still_present_at_warm_up = still_present_at_warm_up[
(still_present_at_warm_up["arrival_time"] < warm_up)
& (
still_present_at_warm_up["depart_time"].isna()
| (still_present_at_warm_up["depart_time"] >= warm_up)
)
]
still_present_at_warm_up| patient | arrival_time | depart_time | |
|---|---|---|---|
| 8 | 9 | 26.653863 | 113.808649 |
| 9 | 10 | 40.737793 | 117.867247 |
| 10 | 11 | 71.026558 | 132.225321 |
| 11 | 12 | 87.458700 | 136.038252 |
| 12 | 13 | 87.465138 | 143.522826 |
| 13 | 14 | 98.810612 | 156.491653 |
| 14 | 15 | 99.173100 | 166.212785 |
These are the patients genuinely still in the system - mid-queue or mid-treatment - at the moment warm_up ends. Any correct animation of the post-warm-up period has to show them.
The problem: filtering the log breaks the animation
The obvious way to discard a warm-up period is event_log[event_log["time"] >= warm_up]. reshape_for_animations warns immediately if you try it:
naive_filtered_log = event_log[event_log["time"] >= warm_up]
naive_reshaped = reshape_for_animations(
naive_filtered_log,
every_x_time_units=10,
limit_duration=g.sim_duration,
entity_col_name="patient",
run_col_name="run",
)/tmp/ipykernel_4853/3218288777.py:3: UserWarning: 7 entities (9, 13, 10, 14, 11, ...) have events in the event log but no 'arrival' event, so they will be missing from every frame of the animation.
vidigi works out who is present at each snapshot from the arrival and departure rows, so an entity without an arrival is never drawn.
The usual cause is discarding a warm-up period by filtering the log, e.g. `event_log[event_log['time'] >= warm_up]`, which removes the arrival rows of everyone already in the system - including entities that are still queuing.
To skip a warm-up period, pass the whole event log and set `warm_up` to the end of the warm-up instead. That trims the animation window without discarding the history it needs.
naive_reshaped = reshape_for_animations(The warning names exactly why: filtering removed the arrival row of every patient who arrived before warm_up but is still in the system - reshape_for_animations works out who is present at each snapshot from the arrival and departure rows, so an entity with no arrival row is never drawn, in any frame, not just the first one:
present_anywhere = set(naive_reshaped["patient"].dropna().unique())
missing = set(still_present_at_warm_up["patient"]) - present_anywhere
f"{len(missing)} of {len(still_present_at_warm_up)} genuinely-present patients never appear anywhere in the naively-filtered animation"'7 of 7 genuinely-present patients never appear anywhere in the naively-filtered animation'The fix: warm_up=
Pass the whole event log instead, and let warm_up= trim the animation window - the history before it is kept around just long enough to work out who is present when the window opens:
proper_reshaped = reshape_for_animations(
event_log,
every_x_time_units=10,
limit_duration=g.sim_duration,
entity_col_name="patient",
run_col_name="run",
warm_up=warm_up,
)
first_frame = proper_reshaped[proper_reshaped["snapshot_time"] == warm_up]
sorted(first_frame["patient"].dropna().unique())[np.int64(9),
np.int64(10),
np.int64(11),
np.int64(12),
np.int64(13),
np.int64(14),
np.int64(15)]Every one of the patients identified above is there, in the very first frame - no warning, nothing missing.
animate_activity_log takes the same warm_up= and forwards it straight through, so the fix is identical whichever level you’re working at:
event_position_df = create_event_position_df(
[
EventPosition(event="arrival", x=50, y=300, label="Arrival"),
EventPosition(
event="treatment_wait_begins", x=205, y=275, label="Waiting for Treatment"
),
EventPosition(
event="treatment_begins",
x=205,
y=175,
resource="n_cubicles",
label="Being Treated",
),
EventPosition(event="depart", x=270, y=70, label="Exit"),
]
)
fig = animate_activity_log(
event_log=event_log,
event_position_df=event_position_df,
entity_col_name="patient",
scenario=g(),
setup_mode=False,
every_x_time_units=10,
warm_up=warm_up,
include_play_button=True,
resource_icon_size=15,
text_size=20,
entity_icon_size=13,
gap_between_entities=6,
gap_between_queue_rows=25,
gap_between_resource_rows=25,
plotly_height=600,
frame_duration=200,
plotly_width=1000,
override_x_max=300,
override_y_max=500,
limit_duration=g.sim_duration,
wrap_queues_at=25,
step_snapshot_max=125,
time_display_units="dhm",
display_stage_labels=False,
add_background_image="https://raw.githubusercontent.com/Bergam0t/vidigi/refs/heads/main/examples/example_1_simplest_case/Simplest%20Model%20Background%20Image%20-%20Horizontal%20Layout.drawio.png",
)
figThe animation opens already populated - a queue and a full set of cubicles, not the empty system the model actually started in - and runs from t=100 to the end of the run.
snapshot_alignment: choosing where the grid starts
warm_up=100 and every_x_time_units=10 divide evenly, so there’s only one sensible snapshot grid. That’s not always true - with every_x_time_units=15, 100 isn’t a multiple of 15, and snapshot_alignment decides what happens:
from_warm_up = reshape_for_animations(
event_log,
every_x_time_units=15,
limit_duration=g.sim_duration,
entity_col_name="patient",
run_col_name="run",
warm_up=warm_up,
snapshot_alignment="warm_up",
)
from_run_start = reshape_for_animations(
event_log,
every_x_time_units=15,
limit_duration=g.sim_duration,
entity_col_name="patient",
run_col_name="run",
warm_up=warm_up,
snapshot_alignment="run_start",
)
(
sorted(from_warm_up["snapshot_time"].unique())[:6],
sorted(from_run_start["snapshot_time"].unique())[:6],
)([np.int64(100),
np.int64(115),
np.int64(130),
np.int64(145),
np.int64(160),
np.int64(175)],
[np.int64(105),
np.int64(120),
np.int64(135),
np.int64(150),
np.int64(165),
np.int64(180)])snapshot_alignment="warm_up" (the default) puts the first frame exactly on the boundary, 100, so the animation opens showing the system precisely as the warm-up ends. snapshot_alignment="run_start" keeps the grid that would have applied with no warm-up at all - multiples of 15 from 0 - and simply drops the frames before 100, so the first surviving one is 105. Frame times then stay the same round numbers a no-warm-up run would have used, at the cost of the first frame no longer landing exactly on the cutoff. The two are identical whenever warm_up happens to be a multiple of every_x_time_units, as it was above.
Closing note
This and feat_warm_up.ipynb answer two different questions about the same idea. That notebook is about choosing a warm-up length for reported statistics, using plot_warm_up_diagnostic/welch_moving_average on a duration or queue-length series. This one assumes a length has already been chosen, and is about applying it to an animation without silently losing entities in the process - a mistake the new warning above exists specifically to catch.
examples/v2_release_additions.ipynb tours the rest of what shipped alongside this in 2.0.0.