import random
import plotly.io as pio
import simpy
from sim_tools.distributions import Exponential, Lognormal
from vidigi import analysis
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: A Tour of vidigi 2.0.0’s Analytics & Plotting Additions
Before 2.0.0, vidigi’s non-animation output was four methods: plot_entity_timeline, generate_dfg, plot_metric_bar (bare bars, no uncertainty), and plot_queue_size. 2.0.0 adds a proper analysis surface on top of the same event logs: duration distributions, resource utilisation, confidence intervals across replications, a warm-up diagnostic, a replication-count diagnostic, a metric-vs-arrival-time diagnostic, an outlier-run diagnostic, a rare-event rate estimate, a two-scenario comparison highlighter, and a per-replication metric boxplot with configurable threshold shading - plus a few smaller quality-of-life additions to the logging and animation APIs themselves.
Every vidigi.plots function is a thin wrapper over a matching vidigi.analysis function that returns the numbers alone, for tables and reports; TrialLogger gets one-line delegating methods for both, which is the route used throughout this notebook. Three of these features - warm-up diagnostics, replication-count diagnostics, and metric-vs-arrival-time - already have their own dedicated notebook, so they get a short, real demonstration here with a link to the full treatment rather than being re-explained from scratch. The animation picked up new arguments too: feat_animation_warm_up.ipynb covers discarding a warm-up period from an animation, and the queue_direction section near the end here is the short version of feat_queue_direction.ipynb.
Model setup
The same single-resource clinic model used by feat_warm_up.ipynb, feat_replication_analysis.ipynb and feat_metric_vs_arrival_time.ipynb - patients arrive, wait for one of 4 treatment cubicles, are treated, and leave. Reusing it means every number below is directly comparable with those three notebooks (the ~78% cubicle utilisation quoted there is the same figure this notebook derives independently in the resource-utilisation section).
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 = 20class Patient:
"""Class defining details for a patient entity"""
def __init__(self, p_id):
self.id = p_idclass Model:
def __init__(self, run_number, n_cubicles=None):
self.env = simpy.Environment()
self.run_number = run_number
# n_cubicles= is an override for the scenario-comparison section further
# down; every other cell in this notebook omits it, so self.n_cubicles is
# a verified no-op equal to g.n_cubicles.
self.n_cubicles = n_cubicles if n_cubicles is not None else g.n_cubicles
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=self.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)class Trial:
def __init__(self, n_cubicles=None):
self.n_cubicles = n_cubicles
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, n_cubicles=self.n_cubicles)
my_model.run()
self.all_event_logs.append(my_model.logger)clinic_trial = Trial()
trial_logs = TrialLogger(
clinic_trial.all_event_logs, scenario=g(), label="base scenario"
)
event_log = trial_logs.to_dataframe()
trial_logs.summary(){'number_of_runs': 20, 'label': 'base scenario', 'scenario_attached': True}Getting the raw numbers: event_durations / get_event_durations
Everything else in this notebook - distributions, bar charts, warm-up and replication diagnostics, metric-vs-arrival-time - is built on one function: event_durations, which pairs two events per entity and returns a duration, one row per pairing, with no aggregation. TrialLogger.get_event_durations is the same thing called on the trial’s combined log. It replaces an older pivot-based calculation that raised on any entity revisiting a step (a rework loop); this one supports that via match="first"|"last"|"occurrence".
wait_durations = trial_logs.get_event_durations(
"treatment_wait_begins", "treatment_begins"
)
wait_durations.head()| entity_id | run_number | pathway | occurrence | first_time | second_time | duration | |
|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | <NA> | 0 | 0.000000 | 0.000000 | 0.0 |
| 1 | 2 | 1 | <NA> | 0 | 19.233669 | 19.233669 | 0.0 |
| 2 | 3 | 1 | <NA> | 0 | 37.923186 | 37.923186 | 0.0 |
| 3 | 4 | 1 | <NA> | 0 | 57.001274 | 57.001274 | 0.0 |
| 4 | 5 | 1 | <NA> | 0 | 59.239628 | 59.239628 | 0.0 |
wait_durations["duration"].describe()count 7475.000000
mean 8.323483
std 12.148811
min 0.000000
25% 0.000000
50% 1.686436
75% 13.231493
max 76.838840
Name: duration, dtype: float64entity_id, run_number, pathway and occurrence come along for free, and a pairing with no matching second event (still queuing when the run ends) is kept with duration = NaN rather than dropped, by default (keep_incomplete=True).
One summary number, pooled or per-replication: get_event_duration_stat / get_event_duration_ci
get_event_durations gives every duration; often you want just one number. get_event_duration_stat reduces them to a single statistic. Its default across="entities" pools every entity’s duration across every replication - one big sample, run boundaries ignored - exactly as every prior release. across="runs" (new in 2.0.0) instead computes the statistic within each run and averages those, weighting every replication equally rather than by how many patients it happened to see. get_event_duration_ci (also new) goes one step further, returning that per-replication mean with a Student’s t confidence interval across the 20 runs - the headline “what is this number, and how sure are we” figure, and the numbers-only twin of the plot_metric_bar(across="runs", error_bars="ci") bar in the next section.
See Choosing how to summarise across replications for which one a given question wants, and feat_trial_logger.ipynb for the full treatment.
pooled = trial_logs.get_event_duration_stat("treatment_wait_begins", "treatment_begins")
per_run = trial_logs.get_event_duration_stat(
"treatment_wait_begins", "treatment_begins", across="runs"
)
ci = trial_logs.get_event_duration_ci("treatment_wait_begins", "treatment_begins")
print(f"pooled over all patients : {pooled}")
print(f"mean of per-run means : {per_run}")
print(
f"per-run mean with 95% CI : {ci.mean:.2f} ({ci.lower:.2f} to {ci.upper:.2f}, n={ci.n})"
)pooled over all patients : 8.32
mean of per-run means : 8.24
per-run mean with 95% CI : 8.24 (6.80 to 9.68, n=20)Seeing the shape of a duration: plot_duration_distribution
A single statistic hides the shape - two systems with the same mean wait can look very different once you see the spread. plot_duration_distribution draws it as a histogram, box, violin, ECDF, or (given split_by) a ridgeline or heatmap comparing many groups at once. All six kind= options, and split_by="run"|"pathway", are demonstrated in full in feat_trial_logger.ipynb - here’s the one most useful for this model, a violin per run:
trial_logs.plot_duration_distribution(
"treatment_wait_begins",
"treatment_begins",
kind="violin",
split_by="run",
title="Treatment wait, by run",
)Uncertainty on a bar chart: plot_metric’s across=/error_bars=
plot_metric (new in 2.0.0, kind="bar" by default) draws a bar chart with a sense of how much a statistic varies between replications, rather than a bare number. across="runs" computes the statistic separately within each run first, then error_bars="ci" draws a confidence interval over those per-run values (the only statistically valid unit to interval over - entities within a run are correlated, replications are not). show_runs=True overlays each run’s own value as a point. This replaces plot_metric_bar, now deprecated (it remains available unchanged until 3.0, but emits a DeprecationWarning) - plot_metric is built on plotly.graph_objects rather than plotly.express, which is what lets it also offer kind="box"/"violin" and highlight_bands below without changing what plot_metric_bar’s own **kwargs means to existing callers. Full coverage, including every error_bars option, in feat_trial_logger.ipynb:
fig = trial_logs.plot_metric(
[
{
"first_event": "treatment_wait_begins",
"second_event": "treatment_begins",
"label": "Treatment wait",
}
],
what="mean",
across="runs",
error_bars="ci",
show_runs=True,
)
fig.update_layout(title="Mean treatment wait, with 95% CI across runs", width=700)
fig.show()Boxplots of per-replication metrics: plot_metric’s kind= and highlight_bands
The bar above collapses each run down to one confidence interval. plot_metric’s kind="box"/"violin" draws the full distribution of per-run values instead - this is the “boxplots of metrics” gap from #153 that plot_duration_distribution’s own kind="box" doesn’t close, since that one only ever plots raw per-entity durations (optionally split by run/pathway), never a box built from a per-run statistic - it requires across="runs", the same as error_bars above.
highlight_bands draws shaded threshold zones behind the chart - useful for a target/acceptable range, or an escalation threshold, read directly off the plot rather than eyeballed against an axis. Each band takes lower/upper (either optional - a missing bound extends to the plotted data’s range), a colour, and an optional label for a legend entry:
fig = trial_logs.plot_metric(
[
{
"first_event": "treatment_wait_begins",
"second_event": "treatment_begins",
"label": "Treatment wait",
}
],
kind="box",
across="runs",
show_runs=True,
highlight_bands=[
{"upper": 6, "colour": "green", "label": "target (<6 min)"},
{"lower": 12, "colour": "red", "label": "concern (>=12 min)"},
],
)
fig.update_layout(title="Per-run treatment wait against target/concern thresholds", width=700)
fig.show()5 of the 20 runs land in the green “target” zone (under 6 minutes), 2 cross into the red “concern” zone (12 minutes or more, including run 10’s 15.96 - the same run that came closest to the outlier fence, without crossing it, in get_outlier_runs above), and the remaining 13 sit in between. The box’s own quartiles and the individual run points (show_runs=True) are still there underneath the bands - highlight_bands adds a reference, it doesn’t replace the distribution itself.
How busy was this resource?
Before 2.0.0 there was no way to answer this at all - the events were there (resource_use/resource_use_end pairs on resource_id), but nothing computed busy time or utilisation from them. Three new vidigi.analysis functions build on top of resource_use_intervals, which pairs the raw resource_use/resource_use_end rows into one interval per bout of use:
intervals = analysis.resource_use_intervals(event_log)
intervals.head()/tmp/ipykernel_5572/2955262641.py:1: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = analysis.resource_use_intervals(event_log)| run_number | entity_id | event | resource_id | start | end | censored | busy_time | |
|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | treatment_begins | 1.0 | 0.000000 | 26.039440 | False | 26.039440 |
| 1 | 1 | 2 | treatment_begins | 2.0 | 19.233669 | 39.185177 | False | 19.951508 |
| 2 | 1 | 3 | treatment_begins | 3.0 | 37.923186 | 66.365742 | False | 28.442556 |
| 3 | 1 | 4 | treatment_begins | 4.0 | 57.001274 | 86.535121 | False | 29.533847 |
| 4 | 1 | 5 | treatment_begins | 1.0 | 59.239628 | 75.897422 | False | 16.657793 |
An entity still holding a resource when the trial’s window ends is censored by default (unclosed="censor"), not dropped - its interval is clipped to the window end rather than discarded, since dropping it would understate utilisation exactly when it matters most (a system that is still busy at the end of a run is disproportionately a congested one). Here that affects a small minority of bouts:
print(
f"{intervals['censored'].sum()} of {len(intervals)} bouts were still open when their run ended"
)66 of 7475 bouts were still open when their run endedresource_occupancy_over_time is the resource equivalent of plot_queue_size’s data - how many units were busy at regular snapshots, computed exactly via a +1/-1 sweep over the intervals above rather than a per-snapshot scan:
occupancy = analysis.resource_occupancy_over_time(event_log, every_x_time_units=50)
occupancy.head()/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:2165: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(| run_number | event | snapshot_time | count | |
|---|---|---|---|---|
| 0 | 1 | treatment_begins | 0 | 1 |
| 1 | 1 | treatment_begins | 50 | 1 |
| 2 | 1 | treatment_begins | 100 | 1 |
| 3 | 1 | treatment_begins | 150 | 3 |
| 4 | 1 | treatment_begins | 200 | 4 |
resource_utilisation: four ways to supply capacity
Busy time and how many units were in use on average (mean_in_use) need no capacity at all. utilisation (mean_in_use / capacity) does, and there are four ways to supply it, in precedence order. All four resolve to the same answer here, since they’re describing the one real model:
- A - an explicit dict, the simplest route when you already know the numbers:
ru_a = analysis.resource_utilisation(
event_log, resource_capacities={"treatment_begins": 4}
)
ru_a["utilisation"].agg(["mean", "min", "max"])/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(mean 0.775829
min 0.726991
max 0.834697
Name: utilisation, dtype: float64- B - a
scenarioplus a step-to-attribute mapping, so capacity is read directly off the same scenario the model itself uses (no number to keep in sync by hand).scenariocan be that parameter object or a plain dict (scenario={"n_cubicles": 4}) - the mapping’s values are resolved as attributes or dict keys interchangeably:
resource_map = {"treatment_begins": "n_cubicles"}
ru_b = analysis.resource_utilisation(event_log, scenario=g(), resource_map=resource_map)
ru_b["utilisation"].agg(["mean", "min", "max"])/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(mean 0.775829
min 0.726991
max 0.834697
Name: utilisation, dtype: float64- C - a
scenarioplus anevent_position_df, reusing the sameresource=column the animation functions already read capacity from - handy when one already exists for the animation and you don’t want to write a second mapping. Pass such anevent_position_dfto an animation without ascenarioand vidigi now warns rather than silently drawing no resource icons:
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,
label="Being Treated",
resource="n_cubicles",
),
EventPosition(event="depart", x=270, y=70, label="Exit"),
]
)
ru_c = analysis.resource_utilisation(
event_log, scenario=g(), event_position_df=event_position_df
)
ru_c["utilisation"].agg(["mean", "min", "max"])/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(mean 0.775829
min 0.726991
max 0.834697
Name: utilisation, dtype: float64- D -
capacity="infer", for a log with no scenario at all (e.g. fromciwor a CSV): capacity is estimated as the number of distinctresource_ids seen for a step. It’s always a lower bound - a unit that was never used is invisible - and always warns. It happens to recover the true capacity exactly here, because every one of the 4 cubicles gets used at some point across 3000 time units:
ru_d = analysis.resource_utilisation(event_log, capacity="infer")
ru_d["capacity"].unique()/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(
/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1965: UserWarning: capacity='infer' estimates each step's capacity as the number of distinct resource_id values seen for it in the log. This is a *lower bound* even in the best case - a resource unit that was never used is invisible, so inferred utilisation is biased upward for underloaded resources. It can be far worse than a small undercount: if resource_id does not identify individual physical units (e.g. every grant of a plain simpy.Resource pool logs the same id, or none at all), this infers a capacity of 1 regardless of the real capacity, and utilisation is inflated by roughly that factor. Pass resource_capacities=, or scenario= with resource_map= or event_position_df=, for an exact answer.
capacities = _resolve_resource_capacities(array([4])Across 20 replications, mean cubicle utilisation is a fairly stable ~78% (bouncing between roughly 73% and 84% run-to-run) - the same figure feat_warm_up.ipynb quotes for this model, derived independently here rather than carried over.
trial_logs was built with scenario=g() attached (see Carrying the scenario with the trial below), so the TrialLogger resource-utilisation methods from here on read capacity straight from it - resource_map is all they need, no scenario= on each call. The vidigi.analysis free functions above still take scenario= explicitly; only the TrialLogger methods have somewhere to read it from.
by="resource" breaks utilisation down per physical unit instead of per step - useful for spotting one persistently busier cubicle. It needs resource_id to be unique across the whole log, not just within one step; resource_col_name=None on the TrialLogger methods (used below) auto-detects unique_resource_id when it’s present, which is what VidigiStore(..., label=...) logs alongside the plain resource_id used above - the collision-safe route, covered in full in feat_trial_logger.ipynb:
by_resource = trial_logs.get_resource_utilisation(by="resource", resource_col_name=None)
by_resource.groupby("resource_id")["utilisation"].mean()/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(resource_id
treatment_cubicle_1 0.779718
treatment_cubicle_2 0.777525
treatment_cubicle_3 0.777520
treatment_cubicle_4 0.768555
Name: utilisation, dtype: float64The two charts
plot_resource_utilisation is the bar-chart counterpart to resource_utilisation (mean across runs, CI error bar, error_bars="ci"/show_runs=True by default since this is new rather than extracted from older bar-only behaviour), with a dashed reference line at 100% - utilisation can never legitimately exceed it, so a bar crossing it is diagnostic of a capacity or logging problem. plot_resource_utilisation_over_time is the resource equivalent of plot_queue_size, one step function per run plus a bold mean. highlight_bands (new in 2.0.0) works on both charts, and spans every facet when plot_resource_utilisation_over_time plots more than one step - see feat_trial_logger.ipynb for a multi-step example:
fig = trial_logs.plot_resource_utilisation(resource_map=resource_map)
fig.update_layout(title="Treatment cubicle utilisation")
fig.show()/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(fig = trial_logs.plot_resource_utilisation_over_time(
every_x_time_units=25,
as_proportion=True,
resource_map=resource_map,
highlight_bands=[{"lower": 0.9, "colour": "red", "label": "near capacity (>=90%)"}],
)
fig.update_layout(
width=900,
height=500,
title="Treatment cubicle occupancy over time (as a proportion of capacity)",
)
fig.show()/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:2165: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (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: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(Occupancy at each step: activity_occupancy_stats
plot_queue_size and resource_occupancy_over_time both answer how many entities were at this step over time. activity_occupancy_stats collapses that to one row per step - the mean, minimum, maximum and median number present - for every queue step and every resource step at once. across_runs="average" (below) reports the figure expected per replication; across_runs="pool" takes one statistic over every run and snapshot together, so the maximum becomes the worst queue seen in any run.
It is the slow one in this notebook - the queue figures are rebuilt the same way animate_activity_log reconstructs each frame, once per run - so every_x_time_units is worth turning up on a long trial.
analysis.activity_occupancy_stats(event_log, every_x_time_units=50)/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:2165: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 3000 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(| event | kind | mean_occupancy | min_occupancy | max_occupancy | median_occupancy | |
|---|---|---|---|---|---|---|
| 0 | treatment_wait_begins | queue | 1.060656 | 0.0 | 7.25 | 0.2 |
| 1 | treatment_begins | resource | 3.032787 | 0.0 | 4.00 | 3.5 |
This table is what EventLogger.generate_dfg(occupancy_metrics=True) merges onto a process map’s nodes, so a directly-follows graph can show where the queue built up and how close each resource ran to capacity - see feat_process_maps.ipynb.
Choosing a warm-up length: plot_warm_up_diagnostic
Welch’s procedure, visualised: ensemble-average a series across replications, then smooth it with several window widths, and read off where they agree it’s flattened out. Full derivation and the “trim to 500” recommendation for this exact model are in feat_warm_up.ipynb; here’s the diagnostic itself, on the treatment queue:
fig = trial_logs.plot_warm_up_diagnostic(
series="queue",
event="treatment_wait_begins",
every_x_time_units=25,
windows=(5, 10, 20),
)
fig.update_layout(
width=900, height=500, title="Welch's procedure: treatment queue length"
)
fig.show()Choosing a replication count: plot_replication_analysis
The cumulative mean and its confidence interval as replications accumulate, plus the relative half-width underneath against a dashed 5% reference - the point where it stays below that line is the recommended minimum replication count. feat_replication_analysis.ipynb shows this model needs far more than the 20 replications run here to actually converge - the chart below shows exactly why that’s true at 20 alone:
fig = trial_logs.plot_replication_analysis("treatment_wait_begins", "treatment_begins")
fig.update_layout(width=900, height=650)
fig.show()Spotting a fluke replication: get_outlier_runs
More replications means more chances one of them is a fluke - a run whose mean happens to land well away from the rest, by chance rather than because anything about the model changed. get_outlier_runs flags that directly: it takes the same per-replication values get_replication_precision/plot_replication_analysis are built on, and applies Tukey’s classic 1.5×IQR fence - the same convention error_bars="iqr" on plot_metric_bar/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 (closes #153):
outliers = trial_logs.get_outlier_runs("treatment_wait_begins", "treatment_begins")
outliers| run_number | value | lower_fence | upper_fence | is_outlier | |
|---|---|---|---|---|---|
| 0 | 1 | 9.896811 | -0.53987 | 16.955189 | False |
| 1 | 2 | 6.153969 | -0.53987 | 16.955189 | False |
| 2 | 3 | 10.403102 | -0.53987 | 16.955189 | False |
| 3 | 4 | 12.716393 | -0.53987 | 16.955189 | False |
| 4 | 5 | 4.669498 | -0.53987 | 16.955189 | False |
| 5 | 6 | 10.391688 | -0.53987 | 16.955189 | False |
| 6 | 7 | 6.920161 | -0.53987 | 16.955189 | False |
| 7 | 8 | 5.726887 | -0.53987 | 16.955189 | False |
| 8 | 9 | 6.118741 | -0.53987 | 16.955189 | False |
| 9 | 10 | 15.958124 | -0.53987 | 16.955189 | False |
| 10 | 11 | 5.187044 | -0.53987 | 16.955189 | False |
| 11 | 12 | 5.426378 | -0.53987 | 16.955189 | False |
| 12 | 13 | 6.121907 | -0.53987 | 16.955189 | False |
| 13 | 14 | 9.076684 | -0.53987 | 16.955189 | False |
| 14 | 15 | 8.295589 | -0.53987 | 16.955189 | False |
| 15 | 16 | 6.337659 | -0.53987 | 16.955189 | False |
| 16 | 17 | 11.142540 | -0.53987 | 16.955189 | False |
| 17 | 18 | 11.706501 | -0.53987 | 16.955189 | False |
| 18 | 19 | 4.869758 | -0.53987 | 16.955189 | False |
| 19 | 20 | 7.610382 | -0.53987 | 16.955189 | False |
None of the 20 runs is flagged here - run 10 (15.96) comes closest to the upper fence (16.96) without crossing it, and the rest sit well inside. That’s a perfectly normal, healthy result: not every trial has a fluke replication, and get_outlier_runs returning nothing flagged is the confirmation, not a null result. Widening iqr_multiplier narrows the fence further inward for a stricter check (3.0, Tukey’s “far out” convention, flags only the most extreme cases); there’s no method= for a different fencing rule - matching this package’s other diagnostics, one well-understood default rather than a menu of options to choose between.
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:
fig = trial_logs.plot_outlier_runs("treatment_wait_begins", "treatment_begins")
fig.update_layout(width=900, height=400)
fig.show()The title states the same verdict get_outlier_runs’s table gave above, as a sentence: “0 of 20 runs flagged (Tukey 1.5x IQR fence: -0.54 to 17)” - 16.96 rounded to 3 significant figures. The picture backs it up: every point sits inside the shaded fence bands, packed onto a single row since no two runs are close enough together to need spreading out, and there’s no "outlier run" trace in the legend at all - with nothing flagged, plot_outlier_runs omits that trace entirely rather than drawing an empty one.
A rare-event estimate: event_occurrence_rate
get_event_duration_stat’s "unserved_rate"/"served_rate" (see above) answer a per-entity question: what fraction of patients reached one event but not another. event_occurrence_rate (new in 2.0.0) answers a different, per-run question instead: in what proportion of replications did some event happen at all, at least once - the shape a genuinely rare, binary condition takes (a capacity breach, a specific alarm), not a duration between two events. Its interval is a Wilson score interval, not get_event_duration_ci’s Student’s t - a proportion’s interval needs to stay inside [0, 1] and behave sensibly near 0 or 1, exactly where a rare-event rate typically sits (closes #153).
This clinic model has no naturally rare event of its own to demonstrate on, so here’s one derived from a question the resource-utilisation and warm-up sections above didn’t ask: in how many of the 20 replications did the treatment queue ever reach 11 or more patients waiting at once - a genuinely rare “is this a busy day” alarm, not just typical variation? queue_size_over_time (used inside plot_warm_up_diagnostic earlier, called directly here) gives the snapshot-by-snapshot queue length needed to answer it:
import pandas as pd
queue_lengths = analysis.queue_size_over_time(
event_log, ["treatment_wait_begins"], limit_duration=g.sim_duration,
every_x_time_units=5,
)
per_run_max = queue_lengths.groupby("run_number")["count"].max()
overflow_runs = sorted(per_run_max[per_run_max >= 11].index)
print(f"queue reached 11+ at some snapshot in runs: {overflow_runs}")
queue_overflow_log = pd.DataFrame({"run_number": overflow_runs, "event": "queue_overflow"})
rate = analysis.event_occurrence_rate(
queue_overflow_log, "queue_overflow", n_runs=g.number_of_runs
)
ratequeue reached 11+ at some snapshot in runs: [3, 10]ProportionEstimate(proportion=0.1, lower=np.float64(0.027866481213768224), upper=np.float64(0.3010336452284873), n_runs=20, n_occurred=2, ci_level=0.95, method='wilson')A queue of 11+ waiting patients at once happened in only 2 of the 20 runs (10%) - genuinely rare. The Wilson interval is wide (2.8% to 30.1%) despite that clean 10% point estimate, because it’s built from only 20 observations - a reminder that a rare-event rate estimated from a modest number of replications is itself imprecise, the same “more replications needed” message plot_replication_analysis gives for a duration statistic, just for a proportion instead. TrialLogger.get_event_occurrence_rate(event_name) is the equivalent one-line delegator for a model that logs the rare condition directly (e.g. via log_custom_event the moment a real capacity-breach alarm fires) rather than one derived after the fact like this.
Does it matter when you arrived? plot_metric_vs_arrival_time
A third question, distinct from both of the above: not “how many replications” or “how much warm-up”, but whether a metric drifts within one run depending on when the entity that produced it arrived - a non-stationary arrival process, or a time-of-day load effect. Full treatment, including rolling_window vs rolling_time smoothing and colour_by, is in feat_metric_vs_arrival_time.ipynb. highlight_bands (new in 2.0.0, same shape as plot_metric’s above) works here too, shading a target zone against arrival time rather than run:
fig = trial_logs.plot_metric_vs_arrival_time(
"treatment_wait_begins",
"treatment_begins",
rolling_time=150,
marker_size=3,
highlight_bands=[{"upper": 5, "colour": "green", "label": "target (<5 min)"}],
)
fig.update_layout(width=900, height=500)
fig.show()Comparing two scenarios: compare_event_duration_stat / plot_event_duration_comparison / compare_resource_utilisation
Everything above describes one scenario. The question a capacity-planning exercise usually asks next is comparative: if the clinic ran with fewer cubicles, would patients notice? compare_event_duration_stat - and its resource-utilisation twin, compare_resource_utilisation - take a second TrialLogger and answer exactly that: the scenario comparison highlighter from #153.
Model/Trial above take an optional n_cubicles= override (a verified no-op when omitted - every run so far used the default g.n_cubicles, 4). A second trial with 3 cubicles instead of 4 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 inter-arrival and treatment-time streams are seeded from run_number * g.random_number_set, independent of n_cubicles - so this is a matched (common-random-numbers) comparison, the standard way to keep the variance of a difference down when comparing two scenarios’ replications pairwise. The reduced scenario gets its own small subclass of g, rather than mutating g.n_cubicles directly, so trial_logs.scenario.n_cubicles (attached to the first trial) keeps reading 4:
class GFewerCubicles(g):
n_cubicles = 3
reduced_trial = Trial(n_cubicles=GFewerCubicles.n_cubicles)
reduced_trial_logs = TrialLogger(
reduced_trial.all_event_logs, scenario=GFewerCubicles(), label="3 cubicles"
)
reduced_trial_logs.summary(){'number_of_runs': 20, 'label': '3 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, "treatment_wait_begins", "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: 8.2
3 cubicles mean wait: 108.6
delta: +100.4 (+1219%)
CIs overlap: False Welch's t-test p-value: 6.401e-09plot_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 - useful for seeing at a glance which scenario’s bar actually lands inside an acceptable range:
fig = trial_logs.plot_event_duration_comparison(
reduced_trial_logs,
"treatment_wait_begins",
"treatment_begins",
highlight_bands=[{"upper": 15, "colour": "green", "label": "acceptable (<15 min)"}],
)
fig.update_layout(width=700, height=500)
fig.show()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 compare_replication_values to compare one specific step or resource instead. Each trial resolves capacity from its own attached scenario (g() for trial_logs, GFewerCubicles() for reduced_trial_logs), so only resource_map needs passing here, exactly as in the single-scenario resource section above:
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}"
)base scenario utilisation: 78%
3 cubicles utilisation: 97%
CIs overlap: False Welch's t-test p-value: 4.748e-22/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(
/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 60 resource_use row(s) were still open at the end of the window - censored at 2999.936595106806 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(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)
fig.show()/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(
/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 60 resource_use row(s) were still open at the end of the window - censored at 2999.936595106806 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(Same verdict as the numbers above, now a title: “95% CIs do not overlap (p=4.75e-22, Welch’s t-test)”, with two bars - 78% for base scenario, 97% for 3 cubicles - and error bars visibly clear of each other. It’s the same chart mechanics as plot_event_duration_comparison above, just fed by compare_resource_utilisation instead of compare_event_duration_stat.
Both comparisons land the same way: three cubicles is unambiguously worse. Mean treatment wait balloons from 8.2 to 108.6 time units (+1219%) and cubicle utilisation climbs from 78% to 97% - both with non-overlapping 95% CIs and a p-value indistinguishable from zero. None of that is a surprise for a queueing system this close to saturation (97% utilisation leaves almost no slack to absorb variation in arrivals or treatment times), but it’s exactly the kind of question compare_event_duration_stat/compare_resource_utilisation are for: turning “would patients notice?” into a number, a confidence interval, and a plain-language verdict, rather than eyeballing two separate plot_metric_bar charts and guessing.
Keeping an entity timeline: plot_entity_timeline’s new return_fig=
EventLogger.plot_entity_timeline has always shown one entity’s own journey through the model - useful for debugging a specific case. Until 2.0.0 it only ever called fig.show() and returned None; there was no way to keep the figure to restyle or export it. return_fig=True returns it instead:
fig = clinic_trial.all_event_logs[0].plot_entity_timeline(5, return_fig=True)
fig.update_layout(title="Entity 5's journey through the clinic (run 1)")
fig.show()
# fig.write_image("entity_5_timeline.png") # now possible, since fig is a real FigureThe default is still False - existing scripts that rely on plot_entity_timeline displaying itself keep working unchanged. The default is planned to flip to True at vidigi 3.0.
Animating from a logger directly: logger.animate_activity_log() and run_number
The animation sections that follow build a one-run DataFrame with clinic_trial.all_event_logs[0].to_dataframe() before calling the animate_activity_log function. As of 2.0.0 that conversion is optional - an EventLogger or TrialLogger can be animated directly, either passed to the function as event_log= or, new, by calling animate_activity_log() (or reshape_for_animations()) as a method on the logger itself - no import, no .to_dataframe():
# an EventLogger animates its own single run - event_position_df is the only argument it needs
clinic_trial.all_event_logs[0].animate_activity_log(event_position_df, ...)For a TrialLogger, run_number= picks the replication, and the scenario attached at construction (here g(), see Carrying the scenario with the trial below) is reused - so the call below passes neither a DataFrame nor scenario=, and resource-availability icons still appear:
# No import, no .to_dataframe(), no scenario= - the TrialLogger has all three
trial_logs.animate_activity_log(
event_position_df,
run_number=1,
every_x_time_units=50,
limit_duration=1000,
plotly_height=450,
plotly_width=1000,
)run_number applies only to a TrialLogger - passing a multi-run TrialLogger without it, or run_number alongside a DataFrame or an EventLogger, raises a ValueError. The animate_activity_log(event_log=...) function form accepts all three inputs too; a plain DataFrame works exactly as before, which is what the animation cells below use so they can also show the function itself.
Carrying the scenario with the trial: scenario= and label=
EventLogger and TrialLogger take an optional scenario= - the parameter object (or plain dict) the runs came from - and a human-readable label=. This trial was built with both, right at the top of the notebook:
trial_logs = TrialLogger(clinic_trial.all_event_logs, scenario=g(), label="base scenario")scenario is the same object-or-dict shape the resource-utilisation helpers accept, so once it is attached get_resource_utilisation, plot_resource_utilisation and plot_resource_utilisation_over_time read capacity straight from it - which is why those calls in the resource section passed only resource_map, never scenario=g(). A per-call scenario= still overrides it. A TrialLogger built from EventLoggers that each carry a scenario / label inherits them (warning if the runs disagree), and vidigi.ciw.event_logger_from_ciw_recs / trial_logger_from_ciw_recs gained matching scenario= / label= arguments.
Both are plain attributes on the logger:
print("label:", trial_logs.label)
print(
"n_cubicles, read back from the attached scenario:", trial_logs.scenario.n_cubicles
)label: base scenario
n_cubicles, read back from the attached scenario: 4Saving a finished trial: to_pickle / read_pickle
Twenty runs of a stochastic model isn’t free to regenerate, and the numbers only mean anything alongside the parameters that produced them. to_pickle() writes the whole trial - every run’s events, the attached scenario, the label - to a single file; TrialLogger.read_pickle() brings it back as a fully working TrialLogger (EventLogger has the same pair).
A logger built with env= (the normal simpy pattern) used to be unpicklable - the live simpy.Environment holds generators. The env is now dropped on pickle, since it is only read while logging, so what you reload is a complete, finished record.
import tempfile
from pathlib import Path
pkl_path = Path(tempfile.mkdtemp()) / "base_scenario_trial.pkl"
trial_logs.to_pickle(pkl_path)
print(f"wrote {pkl_path.name} ({pkl_path.stat().st_size / 1_000:.0f} kB)")
reloaded = TrialLogger.read_pickle(pkl_path)
print("reloaded:", reloaded.summary())
# a real TrialLogger - the scenario came back with it, so utilisation still
# resolves from resource_map alone, matching the figure from before the round trip
reloaded.get_resource_utilisation(resource_map=resource_map)["utilisation"].agg(
["mean", "min", "max"]
)wrote base_scenario_trial.pkl (1625 kB)
reloaded: {'number_of_runs': 20, 'label': 'base scenario', 'scenario_attached': True}/home/runner/work/vidigi/vidigi/src/vidigi/analysis.py:1945: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 2999.8678663039827 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
intervals = resource_use_intervals(mean 0.775829
min 0.726991
max 0.834697
Name: utilisation, dtype: float64Flipping a queue’s build direction: queue_direction
Queues have always built out to the left of their event_position_df anchor, so the front of the line sits at the bottom-right corner. queue_direction="right" - on animate_activity_log, generate_animation and generate_animation_df - mirrors that: the anchor becomes the bottom-left corner and the queue extends rightwards, wrapped rows included. It reads better with entity emojis that face right, and the in-service icons and stage labels move to match. A per-event direction column on event_position_df (or EventPosition(..., direction="right")) overrides it for a single stage.
The default "left" is a verified no-op - existing animations are unchanged. feat_queue_direction.ipynb is the full walkthrough, including the per-event override; here it is on run 1 of the model above, reusing the event_position_df from the resource-utilisation section:
from vidigi.animation import animate_activity_log
single_run_log = clinic_trial.all_event_logs[0].to_dataframe()
animate_activity_log(
event_log=single_run_log,
event_position_df=event_position_df,
scenario=g(),
every_x_time_units=50,
limit_duration=1000,
queue_direction="right",
plotly_height=450,
plotly_width=1000,
)Flipping an icon in place: flip_entity_icons
queue_direction mirrors the layout to suit an icon’s facing direction. Sometimes the layout isn’t free to move - a fixed background image, say - but the icon still faces the wrong way. flip_entity_icons=True mirrors the icon itself instead, leaving every coordinate (and queue_direction) untouched; a per-event flip_icons column on event_position_df (or EventPosition(..., flip_icons=True)) overrides it for a single stage, the same way direction overrides queue_direction.
The default False is a verified no-op. feat_flip_entity_icons.ipynb is the full walkthrough, including the per-event override and getting the CSS to a page that isn’t a live notebook; here it is on the same run and layout as above:
animate_activity_log(
event_log=single_run_log,
event_position_df=event_position_df,
scenario=g(),
every_x_time_units=50,
limit_duration=1000,
flip_entity_icons=True,
plotly_height=450,
plotly_width=1000,
)Giving stage labels room to breathe: stage_label_offset
display_stage_labels=True draws each stage’s label a fixed 10 data units past its EventPosition anchor - the same gap however large entity_icon_size / text_size are set. Turn those up for a bigger, more legible animation and the label starts to visually crowd the icon it’s labelling, because the gap between them never grows to match. stage_label_offset (new in 2.0.0) makes that gap a parameter instead of a hardcoded constant, so it can be widened alongside the icons (#122).
The default 10 is a verified no-op - identical to the previous hardcoded value. Below, entity_icon_size, resource_icon_size and text_size are all turned up to 40; stage_label_offset=30 keeps the labels clear of the larger icons instead of crowding them the way the old fixed gap would:
animate_activity_log(
event_log=single_run_log,
event_position_df=event_position_df,
scenario=g(),
every_x_time_units=50,
limit_duration=1000,
entity_icon_size=40,
resource_icon_size=40,
text_size=40,
stage_label_offset=30,
plotly_height=450,
plotly_width=1000,
)A different cap per step: step_snapshot_max_overrides
step_snapshot_max is a single number for the whole animation, which forces a compromise: set it high enough to show a genuine bottleneck queue honestly and every small queue also renders dozens of icons; set it low and the bottleneck is hidden behind a “+ N more” label the moment it matters. step_snapshot_max_overrides (new in 2.0.0) takes a {event: cap} dict that overrides the scalar for named events only - {"waiting_for_bed": 250} shows one long queue in full while every other step stays capped at step_snapshot_max. Any event not in the dict uses the scalar, so it stays the fallback; a key that matches no event in the log raises a warning, so a misspelt event name isn’t silently ignored.
The default None is a verified no-op, and like step_snapshot_max itself the argument is accepted by reshape_for_animations and generate_animation_df as well as animate_activity_log. The clinic model here has only one queue, so this cell uses a small purpose-built log with two: twelve patients waiting for triage and twelve waiting for a bed. The scalar cap of 5 applies to triage; wait_bed is overridden to 20 and so is drawn in full:
import pandas as pd
two_queue_rows = []
for pid in range(1, 13):
two_queue_rows += [
(0, pid, "arrival_departure", "arrival"),
(0, pid, "queue", "wait_triage"),
(200, pid, "arrival_departure", "depart"),
]
for pid in range(101, 113):
two_queue_rows += [
(0, pid, "arrival_departure", "arrival"),
(0, pid, "queue", "wait_bed"),
(200, pid, "arrival_departure", "depart"),
]
two_queue_log = pd.DataFrame(
two_queue_rows, columns=["time", "entity_id", "event_type", "event"]
)
two_queue_positions = create_event_position_df(
[
EventPosition(event="arrival", x=50, y=300, label="Arrival"),
EventPosition(event="wait_triage", x=380, y=300, label="Waiting for Triage"),
EventPosition(event="wait_bed", x=380, y=180, label="Waiting for a Bed"),
EventPosition(event="depart", x=270, y=70, label="Exit"),
]
)
animate_activity_log(
event_log=two_queue_log,
event_position_df=two_queue_positions,
every_x_time_units=20,
limit_duration=180,
step_snapshot_max=5,
step_snapshot_max_overrides={"wait_bed": 20},
wrap_queues_at=5,
plotly_height=450,
plotly_width=1000,
)Choosing where arrivals come from: spawn_in_from_arrival
A depart row in event_position_df makes every entity move to a set point before it disappears. Arrivals now behave the same way. Previously a brand-new entity flew in from the plot’s top-left corner, because a point that wasn’t in the previous frame has no position for Plotly to move it from; spawn_in_from_arrival (new in 2.0.0, default True) gives it one - the EventPosition(event="arrival", ...) anchor - so a new entity glides in from there instead, the arrival-side mirror of depart. It works by slipping the entity in at the anchor one snapshot early (invisibly, then visibly), so there is a real prior position to animate the move from.
This is a breaking change on the defaults: an animation whose layout gives "arrival" a position - as every example in these docs does - now has its new entities slide in from there. Pass spawn_in_from_arrival=False for the old top-left fly-in. Only entities that arrive at least two snapshots in are affected; anyone already in the system when the animation opens has no earlier frame to enter from and still flies in. It’s independent of step_snapshot_reveal_pop_in above, which does the same job for entities re-emerging from behind a + N more label.
Here the arrival anchor is moved to the top-right corner and labelled “Entrance”, so the effect is unmistakable - patients slide in from the entrance rather than the corner of the plot:
spawn_position_df = event_position_df.copy()
spawn_position_df.loc[spawn_position_df["event"] == "arrival", ["x", "y", "label"]] = [
320,
320,
"Entrance",
]
animate_activity_log(
event_log=single_run_log,
event_position_df=spawn_position_df,
scenario=g(),
every_x_time_units=50,
limit_duration=1000,
spawn_in_from_arrival=True,
plotly_height=450,
plotly_width=1000,
)Beyond emoji: entity_icon_font, entity_colour_by, resource_icon
Every entity icon has, until now, been an emoji - and being colour fonts, emoji ignore textfont.color entirely, so an entity could never be coloured by anything. entity_icon_font switches the icon to an icon font (Font Awesome, Bootstrap Icons, Material Symbols, or any font already on the page) instead - thousands of monochrome glyphs, which entity_colour_by can then colour meaningfully. Below, patients are coloured by the specific treatment cubicle (resource_id) treating them - the same physical units resource_utilisation measured earlier in this notebook - and the cubicles themselves get a small bed icon via resource_icon. That’s a field set per event position (here via the resource_icon column), overriding custom_resource_icon for that one stage, so different resource stages can each carry their own icon; the value is a plain glyph or, as here, an image (bed_icon.svg, drawn via layout.images rather than as text).
entity_icon_font / resource_icon_font need plotly >= 5.23.0 (an icon font resolves a numeric weight, which older plotly rejects); pip install -U plotly if the cell below raises a 'weight' property error. Emoji animations are unaffected.
feat_custom_icons.ipynb is the full derivation, including the four confirmed Plotly bugs found and worked around along the way. That notebook also covers entity_annotation_by: once an icon carries baked-in extra text (a running length-of-stay figure, say), combining that with flip_entity_icons/entity_icon_font needs a second, independently-styled trace rather than more text appended onto the icon itself, since Plotly gives a single icon’s <text> node one font and one transform for the whole node - see the “Annotating an icon with extra text” section of the docs for the full trade-off against appending.
icon_position_df = event_position_df.copy()
icon_position_df.loc[
icon_position_df["event"] == "treatment_begins", "resource_icon"
] = "bed_icon.svg"
animate_activity_log(
event_log=single_run_log,
event_position_df=icon_position_df,
scenario=g(),
every_x_time_units=50,
limit_duration=1000,
entity_icon_font="font-awesome",
custom_entity_icon_list=[""], # fa-walking
entity_colour_by="resource_id",
entity_colour_map={
"1.0": "crimson",
"2.0": "steelblue",
"3.0": "seagreen",
"4.0": "goldenrod",
"nan": "lightgrey", # not currently in treatment
},
gap_between_resources=30,
plotly_height=450,
plotly_width=1000,
)Automatic resource-use logging: VidigiStore(..., logger=...)
The model above logs every resource_use/resource_use_end event by hand, bracketing request() with log_resource_use_start/log_resource_use_end calls. 2.0.0 adds an opt-in logger= constructor parameter to VidigiStore/VidigiPriorityStore that does this automatically: pass entity_id= (and optionally start_event=/end_event=) into request(), and the start event is logged the moment the request is actually granted, the end event the moment the resource is released - no manual logging calls at all. The same model, rewritten:
class AutoLoggingModel:
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.patient_inter_arrival_dist = Exponential(
mean=g.arrival_rate, random_seed=run_number * g.random_number_set
)
self.treat_dist = Lognormal(
mean=g.treat_mean,
stdev=g.treat_var,
random_seed=run_number * g.random_number_set,
)
# logger= is the only change here - VidigiStore now knows where to log to.
self.treatment_cubicles = VidigiStore(
self.env,
num_resources=g.n_cubicles,
label="treatment_cubicle",
logger=self.logger,
)
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")
# entity_id= (and the event names) are all request() needs to auto-log both
# resource_use and resource_use_end - no log_resource_use_start/_end calls.
with self.treatment_cubicles.request(
entity_id=patient.id,
start_event="treatment_begins",
end_event="treatment_complete",
) as req:
yield req
yield self.env.timeout(self.treat_dist.sample())
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)random.seed(1)
auto_model = AutoLoggingModel(1)
auto_model.run()
auto_intervals = analysis.resource_use_intervals(auto_model.logger.to_dataframe())
auto_intervals.head()/tmp/ipykernel_5572/2705768775.py:5: UserWarning: 1 resource_use row(s) were still open at the end of the window - censored at 2984.5873371573853 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
auto_intervals = analysis.resource_use_intervals(auto_model.logger.to_dataframe())| run_number | entity_id | event | resource_id | start | end | censored | busy_time | |
|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | treatment_begins | 1.0 | 0.000000 | 26.039440 | False | 26.039440 |
| 1 | 1 | 2 | treatment_begins | 2.0 | 19.233669 | 39.185177 | False | 19.951508 |
| 2 | 1 | 3 | treatment_begins | 3.0 | 37.923186 | 66.365742 | False | 28.442556 |
| 3 | 1 | 4 | treatment_begins | 4.0 | 57.001274 | 86.535121 | False | 29.533847 |
| 4 | 1 | 5 | treatment_begins | 1.0 | 59.239628 | 75.897422 | False | 16.657793 |
Run 1 uses the same seeds as run 1 of the hand-logged model above, so the two should agree exactly - same start/end times, same resource_ids, for every bout. Comparing them needs one extra thing pinned down first: resource_use_intervals resolves its analysis window’s end from the latest time seen in the log it’s given unless limit_duration= is passed explicitly, and the single run’s own log genuinely ends earlier (its last event) than the 20-run trial’s combined log (the latest event anywhere in the trial) - so the still-open bout at the end of run 1 would otherwise get censored against two different window ends. Pin limit_duration=g.sim_duration on both sides for a fair, exact comparison:
manual_run_1 = analysis.resource_use_intervals(event_log, limit_duration=g.sim_duration)
manual_run_1 = manual_run_1[manual_run_1["run_number"] == 1].reset_index(drop=True)
auto_run_1 = analysis.resource_use_intervals(
auto_model.logger.to_dataframe(), limit_duration=g.sim_duration
).reset_index(drop=True)
manual_run_1.drop(columns="run_number").equals(auto_run_1.drop(columns="run_number"))/tmp/ipykernel_5572/3832256686.py:1: UserWarning: 66 resource_use row(s) were still open at the end of the window - censored at 3000 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
manual_run_1 = analysis.resource_use_intervals(event_log, limit_duration=g.sim_duration)
/tmp/ipykernel_5572/3832256686.py:4: UserWarning: 1 resource_use row(s) were still open at the end of the window - censored at 3000 (unclosed='censor', the default). Pass unclosed='drop' to exclude them instead.
auto_run_1 = analysis.resource_use_intervals(Truelogger= changes nothing about what gets logged, only how much code it takes to log it.
Inspecting a pool like a resource: .count, .capacity, .n_waiting
VidigiStore and VidigiPriorityStore are stores, but a fixed pool built with num_resources= is really being used as a stand-in for simpy.Resource. 2.0.0 adds the read-only properties you would reach for out of that habit:
.count- units currently in use (checked out of the pool).n_waiting- requests queued for a unit (the store analog oflen(simpy.Resource.queue)).num_resources- the pool size: the running total passed tonum_resources=/populate(), later top-uppopulate()calls included
.count is derived from how many units are currently sitting in the store rather than tracked per request, so it stays correct through filter_fn requests, reneging via cancel_get, and VidigiPriorityStore’s direct holder-to-waiter handoff, with no bookkeeping of your own. The invariant 0 <= count <= num_resources is the same one simpy keeps between Resource.count and Resource.capacity.
Everything in this section applies identically to VidigiPriorityStore.
idle_pool = VidigiStore(
simpy.Environment(), num_resources=g.n_cubicles, label="treatment_cubicle"
)
print(f"capacity : {idle_pool.capacity}")
print(f"num_resources : {idle_pool.num_resources}")
print(f"count : {idle_pool.count}")
print(f"n_waiting : {idle_pool.n_waiting}")capacity : 4
num_resources : 4
count : 0
n_waiting : 0.capacity now returns the pool size
Breaking change in 2.0.0 (#87): .capacity used to return the underlying container’s item limit - float("inf") by default, useless as a resource capacity. For a pool built with num_resources= / populate() and no explicit capacity= it now returns the pool size, mirroring simpy.Resource.capacity and equal to .num_resources. A store built with an explicit capacity=, and a bare VidigiStore(env), are unchanged:
print("pooled :", VidigiStore(simpy.Environment(), num_resources=4, label="bed").capacity)
print("bare :", VidigiStore(simpy.Environment()).capacity)
print("explicit:", VidigiStore(simpy.Environment(), num_resources=4, label="bed", capacity=10).capacity)pooled : 4
bare : inf
explicit: 10Watching occupancy during a run
Because .count and .n_waiting read live state, a monitoring process can sample them exactly the way one would sample simpy.Resource.count - the in-simulation equivalent of the resource_occupancy_over_time / plot_queue_size reconstructions earlier in this notebook, which recover the same numbers from the event log after the fact. Here that monitor runs alongside run 1 of the clinic model:
import pandas as pd
monitored = Model(1)
pool_samples = []
def monitor_pool(env, store):
while True:
pool_samples.append(
{"time": env.now, "in_use": store.count, "queued": store.n_waiting}
)
yield env.timeout(50)
monitored.env.process(monitor_pool(monitored.env, monitored.treatment_cubicles))
random.seed(1)
monitored.run()
pool_samples = pd.DataFrame(pool_samples)
pool_samples["in_use"].agg(["mean", "min", "max"])mean 3.066667
min 0.000000
max 4.000000
Name: in_use, dtype: float64Mean cubicle occupancy of ~3.07 of 4, about 77%, is in line with the utilisation section above - derived here live rather than from the log. in_use never exceeds the capacity of 4, the 0 <= count <= num_resources invariant holding throughout the run:
import plotly.graph_objects as go
fig = go.Figure()
fig.add_scatter(
x=pool_samples["time"], y=pool_samples["in_use"],
name="cubicles in use", line_shape="hv",
)
fig.add_scatter(
x=pool_samples["time"], y=pool_samples["queued"],
name="patients queued", line_shape="hv",
)
fig.add_hline(y=g.n_cubicles, line_dash="dash", annotation_text="capacity")
fig.update_layout(
width=900, height=400,
title="Cubicle pool: .count and .n_waiting through run 1",
xaxis_title="time", yaxis_title="patients",
)
fig.show()Returning too many units now raises
Breaking change in 2.0.0: a simpy.Resource never lets you release more units than it has. A VidigiStore pool used to accept an accidental extra put() - a unit returned twice, or one from another pool - and grow silently past num_resources, surfacing only later as a .count RuntimeError. That return now raises ValueError at the call site:
env = simpy.Environment()
tills = VidigiStore(env, num_resources=2, label="till")
held = [tills.get_direct(), tills.get_direct()]
env.run()
for got in held:
tills.put(got.value) # both units back - the pool is now full
try:
tills.put(held[0].value) # the same unit a second time
except ValueError as err:
print(err)VidigiStore is full - all 2 units are already back in the pool, so this return would grow it past num_resources. A unit was probably returned twice, or a unit from another pool was returned here. Pass strict_capacity=False to VidigiStore(...) if the pool is meant to grow.Pass strict_capacity=False to opt out - for a genuinely elastic pool meant to grow via raw put(). .count can then no longer be trusted (the pool size is no longer fixed) and raises its RuntimeError once the store holds more units than were populated:
env = simpy.Environment()
elastic = VidigiStore(env, num_resources=2, label="till", strict_capacity=False)
held = [elastic.get_direct(), elastic.get_direct()]
env.run()
for got in held:
elastic.put(got.value)
elastic.put(held[0].value) # allowed - the pool grows to 3 units
print("units in pool:", len(elastic.items))
try:
elastic.count
except RuntimeError as err:
print(err)units in pool: 3
VidigiStore.count needs the pool size, which is only tracked when units are added via VidigiStore(num_resources=...) or store.populate(...). This store holds more items (3) than were populated that way (2) - populate_store() and hand .put() calls are not counted.The same RuntimeError is why populate_store() - the old free function for seeding a store - is deprecated in 2.0.0 and removed in 3.0: a pool it fills is invisible to .count / .num_resources, and the num_resources= constructor argument and .populate() method now cover the same job. Build the pool with VidigiStore(env, num_resources=N, label=...), as every cell above does, or store.populate(N, label=...) to top an existing one up.
Closing note
Everything above is reachable two ways: as a free function in vidigi.analysis/ vidigi.plots (works on any event log - vidigi.ciw, a CSV, a hand-built frame, not just EventLogger/TrialLogger), or as a one-line delegating method on TrialLogger (used throughout this notebook, since it’s the route most modellers will actually use).
If your model is built with ciw rather than SimPy, that second route is open to you too: vidigi.ciw gained event_logger_from_ciw_recs(recs, node_name_list=...) and trial_logger_from_ciw_recs(list_of_recs, node_name_list=...) in 2.0.0, which build a populated EventLogger / TrialLogger straight from Simulation.get_all_records() output - the same conversion as the long-standing event_log_from_ciw_recs, but returning the logging object instead of a bare DataFrame. Every TrialLogger method demonstrated above then works on a ciw trial unchanged. See the ciw examples for the conversion itself.
For the three diagnostics given only a brief treatment here, the full story is one click away:
- feat_warm_up.ipynb - how much of a run to discard
- feat_replication_analysis.ipynb - how many replications to run
- feat_metric_vs_arrival_time.ipynb - whether a metric drifts within a run
- feat_trial_logger.ipynb - every
plot_duration_distributionkind, everyplot_metric_bar/plot_resource_utilisationerror-bar option, and theresource_idcollision problemlabel=exists to prevent