import random
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
pio.renderers.default = "notebook"Feature Example: Does This Metric Drift Over the Simulation?
Warning
This example has not yet been updated to reflect all of the new features and recommendations from vidigi 2.0.0, which have further simplified the process of adding vidigi to your model and accessing and modifying your animation.
For now, all of the code below will still work, but check out the getting started page for a full guide to the recommended way to use vidigi 2.0.0.
plot_replication_analysis() asks how many replications are enough; plot_warm_up_diagnostic() asks how much warm-up to discard. Neither asks the question this notebook does: within a single run’s steady operation, does a metric’s value depend on when the entity that produced it arrived? A queue that starts empty gives short waits to its earliest arrivals almost by construction - that is the warm-up transient, seen from a different angle - but the same question applies to any non-stationary effect, such as an arrival rate that genuinely changes over the day.
vidigi.analysis.entity_metric_by_arrival() and vidigi.plots.plot_metric_vs_arrival_time() (also TrialLogger.get_entity_metric_by_arrival()/.plot_metric_vs_arrival_time()) answer this directly: a duration plotted against the arrival time of the entity it belongs to, with an optional rolling-average trend line drawn over the top.
This reuses the same single-queue clinic model as feat_warm_up.ipynb and feat_replication_analysis.ipynb, so all three notebooks’ numbers are directly comparable.
Model setup
Show the global parameter class code
class g:
"""
Create a scenario to parameterise the simulation model
Parameters:
-----------
random_number_set: int
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
treat_mean, treat_var: float
Mean and variance of the treatment duration distribution (Lognormal)
arrival_rate: float
Mean of the exponential inter-arrival time distribution
sim_duration: int
The number of time units the simulation will run for
number_of_runs: int
The number of replications
"""
random_number_set = 42
n_cubicles = 4
treat_mean = 25
treat_var = 5
arrival_rate = 8
sim_duration = 3000
number_of_runs = 20Show the patient class code
class Patient:
"""Class defining details for a patient entity"""
def __init__(self, p_id):
self.id = p_idShow the model code
# Identical to feat_warm_up.ipynb's/feat_replication_analysis.ipynb's model - same
# clinic, same parameters, so all three notebooks' numbers are directly comparable.
class Model:
def __init__(self, run_number):
self.env = simpy.Environment()
self.run_number = run_number
self.logger = EventLogger(env=self.env, run_number=self.run_number)
self.patient_counter = 0
self.init_distributions()
self.init_resources()
def init_distributions(self):
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.treat_mean,
stdev=g.treat_var,
random_seed=self.run_number * g.random_number_set,
)
def init_resources(self):
self.treatment_cubicles = VidigiStore(
self.env, num_resources=g.n_cubicles, label="treatment_cubicle"
)
def generator_patient_arrivals(self):
while True:
self.patient_counter += 1
p = Patient(self.patient_counter)
self.env.process(self.attend_clinic(p))
yield self.env.timeout(self.patient_inter_arrival_dist.sample())
def attend_clinic(self, patient):
self.logger.log_arrival(entity_id=patient.id)
self.logger.log_queue(entity_id=patient.id, event="treatment_wait_begins")
with self.treatment_cubicles.request() as req:
treatment_resource = yield req
self.logger.log_resource_use_start(
entity_id=patient.id,
event="treatment_begins",
resource_id=treatment_resource.id,
unique_resource_id=treatment_resource.unique_id,
)
yield self.env.timeout(self.treat_dist.sample())
self.logger.log_resource_use_end(
entity_id=patient.id,
event="treatment_complete",
resource_id=treatment_resource.id,
unique_resource_id=treatment_resource.unique_id,
)
self.logger.log_departure(entity_id=patient.id)
def run(self):
self.env.process(self.generator_patient_arrivals())
self.env.run(until=g.sim_duration)Show the trial class code
class Trial:
def __init__(self):
self.all_event_logs = []
self.run_trial()
def run_trial(self):
for run in range(1, g.number_of_runs + 1):
random.seed(run)
my_model = Model(run)
my_model.run()
self.all_event_logs.append(my_model.logger)clinic_trial = Trial()
trial_logs = TrialLogger(clinic_trial.all_event_logs)
trial_logs.summary(){'number_of_runs': 20, 'label': None, 'scenario_attached': False}Waiting time against arrival time
plot_metric_vs_arrival_time() needs an event pair, just like plot_metric_bar/plot_replication_analysis - here, waiting time from treatment_wait_begins to treatment_begins again. By default it plots against arrival_event="arrival" - deliberately a third, independent event name, not the same as first_event: the question is whether the duration depends on when the entity arrived at the system at all, not on when the specific interval being measured started.
fig = trial_logs.plot_metric_vs_arrival_time(
"treatment_wait_begins",
"treatment_begins",
marker_size=3,
)
fig.update_layout(width=900, height=500)
fig.show()Pooled across all 20 replications this is over 7,000 points - too noisy to read a trend from by eye alone. That is exactly what the rolling-average trend line further down is for; first, the usual source-code and table-behind-the-plot checks.
Checking the claims above against the real implementation
As with the other feature notebooks, the cell below prints the actual entity_metric_by_arrival() source - pulled live from the installed vidigi package via inspect.getsource, not pasted in and liable to drift out of sync with the real code.
Show the entity_metric_by_arrival source, read live from the installed package
import inspect
from vidigi.analysis import entity_metric_by_arrival
print(inspect.getsource(entity_metric_by_arrival))def entity_metric_by_arrival(
event_log: pd.DataFrame,
first_event: str,
second_event: str,
*,
arrival_event: str = "arrival",
match: MatchMode = "first",
entity_col_name: str = "entity_id",
event_col_name: str = "event",
time_col_name: str = "time",
run_col_name: str | None = "auto",
pathway_col_name: str | None = "pathway",
keep_incomplete: bool = True,
) -> pd.DataFrame:
"""
Pair two events per entity, as `event_durations` does, and attach each
entity's arrival time alongside the resulting duration.
Answers a different question from `event_durations` alone: not "how long
did this interval take", but "does this duration vary depending on *when*
the entity arrived at the system" - a non-stationary arrival process or a
time-of-day/load effect, for example. Underlies `plot_metric_vs_arrival_time`.
Parameters
----------
event_log : pandas.DataFrame
Long-format event log, e.g. the output of `EventLogger.to_dataframe()` or
`TrialLogger.to_dataframe()`.
first_event, second_event : str
The two events to pair - see `event_durations`.
arrival_event : str, default="arrival"
The event marking an entity's arrival at the system. Deliberately
independent of `first_event`/`second_event`: it can coincide with
`first_event` (e.g. measuring time from arrival itself), but does not
have to - the arrival time is looked up separately regardless of which
two events the duration itself is measured between.
match : {"first", "last", "occurrence"}, default="first"
How repeated occurrences of `first_event`/`second_event` are paired -
see `event_durations`. Does **not** affect the arrival-time lookup,
which always uses the entity's earliest occurrence of `arrival_event`
regardless of `match` - an entity ordinarily arrives once.
entity_col_name : str, default="entity_id"
Column identifying the entity.
event_col_name : str, default="event"
Column holding the event name.
time_col_name : str, default="time"
Column holding the event time.
run_col_name : str or None, default="auto"
Column identifying which simulation run each row belongs to - see
`event_durations`.
pathway_col_name : str or None, default="pathway"
Column holding the entity's pathway, carried through to the output -
see `event_durations`.
keep_incomplete : bool, default=True
If True (default), rows with no complete `first_event`/`second_event`
pairing are kept, with `NaN` duration - see `event_durations`.
Returns
-------
pandas.DataFrame
One row per matched `(first_event, second_event)` pair - the same
granularity as `event_durations` - with columns ``entity_id``,
``run_number``, ``pathway``, ``occurrence``, ``first_time``,
``second_time``, ``duration``, ``arrival_time``. An entity with a
complete duration pairing but no `arrival_event` recorded in that run
has ``arrival_time = NaN`` rather than being dropped.
Raises
------
ValueError
Everything `event_durations` raises, plus if `arrival_event` is not
present in `event_col_name`.
Notes
-----
The merge joining the duration and arrival frames is keyed on
`(run_number, entity_id)` only, not `occurrence` - every occurrence-row for
one entity (e.g. a rework loop under `match="occurrence"`) shares the same
`arrival_time`, since there is one arrival per entity per run regardless of
how many times the metric's event pair repeats.
See Also
--------
event_durations : The duration pairing this builds on.
vidigi.plots.plot_metric_vs_arrival_time : The matching chart.
Examples
--------
>>> entity_metric_by_arrival(event_log, "treatment_wait_begins", "treatment_begins")
"""
_check_events_present(event_log, event_col_name, arrival_event)
durations = event_durations(
event_log,
first_event,
second_event,
match=match,
entity_col_name=entity_col_name,
event_col_name=event_col_name,
time_col_name=time_col_name,
run_col_name=run_col_name,
pathway_col_name=pathway_col_name,
keep_incomplete=keep_incomplete,
)
arrivals = _first_event_time(
event_log,
arrival_event,
entity_col_name=entity_col_name,
event_col_name=event_col_name,
time_col_name=time_col_name,
run_col_name=run_col_name,
)
merged = durations.merge(arrivals, on=["entity_id", "run_number"], how="left")
return merged[
[
"entity_id",
"run_number",
"pathway",
"occurrence",
"first_time",
"second_time",
"duration",
"arrival_time",
]
].reset_index(drop=True)
The table behind the plot
get_entity_metric_by_arrival() (or the free function vidigi.analysis.entity_metric_by_arrival) returns the same per-entity rows the plot draws, for anyone who wants the table rather than the chart.
table = trial_logs.get_entity_metric_by_arrival(
"treatment_wait_begins", "treatment_begins"
)
table.head()| entity_id | run_number | pathway | occurrence | first_time | second_time | duration | arrival_time | |
|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | <NA> | 0 | 0.000000 | 0.000000 | 0.0 | 0.000000 |
| 1 | 2 | 1 | <NA> | 0 | 19.233669 | 19.233669 | 0.0 | 19.233669 |
| 2 | 3 | 1 | <NA> | 0 | 37.923186 | 37.923186 | 0.0 | 37.923186 |
| 3 | 4 | 1 | <NA> | 0 | 57.001274 | 57.001274 | 0.0 | 57.001274 |
| 4 | 5 | 1 | <NA> | 0 | 59.239628 | 59.239628 | 0.0 | 59.239628 |
first_time is when this specific patient’s wait began - arrival_time is when they arrived at the clinic at all. For this event pair the two happen to be close together (a patient starts waiting almost the moment they arrive), but they are not the same column, and arrival_time is what the chart above plots on its x-axis, not first_time. Pass a different arrival_event= - or a different first_event=/second_event= pair entirely - to ask the same drift question about a metric measured somewhere else in the pathway.
Smoothing the trend: rolling_window vs rolling_time
Two mutually exclusive ways to draw a trend line over the scatter, both symmetric and shrinking - not dropping points - at both edges, so every entity stays visible on the chart:
rolling_window=naverages thennearest points on each side, by arrival order - a count-based window.rolling_time=taverages every point withinttime units on either side of a given arrival time - a genuine time-window, so the number of points contributing varies with how bunched up arrivals happen to be locally, rather than staying fixed.
fig = trial_logs.plot_metric_vs_arrival_time(
"treatment_wait_begins",
"treatment_begins",
rolling_time=150,
marker_size=3,
)
fig.update_layout(width=900, height=500)
fig.show()The trend line tells a clearer story than the raw scatter alone: it climbs from around 4.6 time units right at the start of the run up to roughly 9-10 by about t=400, then spends the rest of the run fluctuating in a noisy band - roughly 5 to 13 - rather than settling onto one flat value (this is the same highly-variable ~78% utilisation queue feat_warm_up.ipynb and feat_replication_analysis.ipynb already found needs a lot of replications to pin down - that noise is real signal about the system, not evidence the transient never ended). The early climb is the warm-up transient, seen here from the arrival-time angle rather than plot_warm_up_diagnostic’s snapshot-grid one - two different diagnostics converging on the same conclusion.
A fixed-count window like rolling_window can look rougher than the time-based version above, even on this fairly steady arrival stream: how much wall-clock time a fixed number of neighbours spans changes with how bunched-up arrivals happen to be locally, whereas rolling_time always spans the same time window. With a mean inter-arrival time of 8, rolling_window=19 covers roughly the same span on average as rolling_time=150 above - the broad climb-then-noisy-plateau shape still comes through, just less smoothly.
fig = trial_logs.plot_metric_vs_arrival_time(
"treatment_wait_begins",
"treatment_begins",
rolling_window=19,
marker_size=3,
)
fig.update_layout(width=900, height=500)
fig.show()colour_by: comparing runs
colour_by="run" (or "pathway") draws one trace per group instead of pooling everything together - matching plot_duration_distribution’s split_by. With all 20 runs at once this would be as unreadable as the very first plot above, so here it is shown against a 3-run subset instead.
first_three_runs = TrialLogger(clinic_trial.all_event_logs[:3])
fig = first_three_runs.plot_metric_vs_arrival_time(
"treatment_wait_begins",
"treatment_begins",
colour_by="run",
marker_size=4,
)
fig.update_layout(width=900, height=500)
fig.show()Each run on its own is noisy enough that its early-vs-late average does not always move in the same direction as the pooled trend line above - one of the three even runs lower later than early, by chance. That is not a contradiction: it is the same lesson feat_replication_analysis.ipynb draws out explicitly - a handful of replications is not enough to pin this system down - showing up again here. The climb is real in the pooled trend precisely because pooling across many runs is what averages that per-run noise away.
warm_up=: excluding the transient
warm_up excludes points by arrival_time - this chart’s x-axis - not by first_time, unlike most other warm_up arguments in vidigi. Filtering by first_time instead could draw “excluded” points to the left of the stated cutoff whenever arrival_event differs from first_event, which would be actively confusing for exactly this diagnostic. Applying feat_warm_up.ipynb’s own choice of warm_up=500 here:
fig = trial_logs.plot_metric_vs_arrival_time(
"treatment_wait_begins",
"treatment_begins",
warm_up=500,
rolling_time=150,
marker_size=3,
)
fig.update_layout(width=900, height=500)
fig.show()The climb is gone: the trend line now starts around 8.2 rather than 4.6, and stays in a flatter (if still noisy) band all the way across - the early, artificially-short waits from the empty-queue start have been excluded, consistent with the warm-up length chosen in the companion notebook.
Closing note
Together with feat_warm_up.ipynb and feat_replication_analysis.ipynb, this covers a third lens on trusting a stochastic simulation’s output: how much of each run to discard, how many runs to run, and now - within one run’s steady operation - whether a metric drifts with when the entity that produced it arrived. All three read off a chart rather than a number the tool decides for you and runs with.