A ciw 2.x example

Warning

Note that this example is written using ciw 2.x

It will not run with 3.x - but could theoretically be adapted to do so

The underlying model code is from Monks, T., Harper, A., & Heather, A. (2023). Towards Sharing Tools, Artefacts, and Reproducible Simulation: a ciw model example (v1.0.1). Zenodo. https://doi.org/10.5281/zenodo.10051494

See here for the adaptation embedded within that repo: https://github.com/Bergam0t/ciw-example-animation/tree/main

In SimPy models, we have to manually add our event logs at various points. However, for Ciw models, we instead can make use of the event_log_from_ciw_recs helper function from vidigi.utils to automatically reshape the logs ciw generates into the correct format for vidigi to work with.

Let’s start by running the model and viewing the logs ciw outputs.

# Import the wrapper objects for model interaction.
from ex_4_ciw_model import Experiment, multiple_replications

from vidigi.ciw import event_logger_from_ciw_recs, trial_logger_from_ciw_recs
from vidigi.utils import EventPosition, create_event_position_df
"""
CiW Implementation of the 111 call centre
Time units of the simulation model are in minutes.
"""
# Imports

import ciw
import numpy as np
import pandas as pd

# Module level variables, constants, and default values

N_OPERATORS = 13
N_NURSES = 9
MEAN_IAT = 100.0 / 60.0

CALL_LOW = 5.0
CALL_MODE = 7.0
CALL_HIGH = 10.0

NURSE_CALL_LOW = 10.0
NURSE_CALL_HIGH = 20.0

CHANCE_CALLBACK = 0.4
RESULTS_COLLECTION_PERIOD = 1000


# Experiment class
class Experiment:
    def __init__(
        self,
        n_operators=N_OPERATORS,
        n_nurses=N_NURSES,
        mean_iat=MEAN_IAT,
        call_low=CALL_LOW,
        call_mode=CALL_MODE,
        call_high=CALL_HIGH,
        chance_callback=CHANCE_CALLBACK,
        nurse_call_low=NURSE_CALL_LOW,
        nurse_call_high=NURSE_CALL_HIGH,
        random_seed=None,
    ):
        self.n_operators = n_operators
        self.n_nurses = n_nurses

        self.arrival_dist = ciw.dists.Exponential(mean_iat)
        self.call_dist = ciw.dists.Triangular(call_low, call_mode, call_high)
        self.nurse_dist = ciw.dists.Uniform(nurse_call_low, nurse_call_high)

        self.chance_callback = chance_callback

        self.init_results_variables()

    def init_results_variables(self):
        self.results = {
            "waiting_times": [],
            "total_call_duration": 0.0,
            "nurse_waiting_times": [],
            "total_nurse_call_duration": 0.0,
        }


# Model code


def get_model(args):
    """
    Build a CiW model using the arguments provided.
    """
    network = ciw.create_network(
        arrival_distributions=[args.arrival_dist, None],
        service_distributions=[args.call_dist, args.nurse_dist],
        routing=[[0.0, args.chance_callback], [0.0, 0.0]],
        number_of_servers=[args.n_operators, args.n_nurses],
    )
    return network


# Model wrapper functions


def single_run(experiment, rc_period=RESULTS_COLLECTION_PERIOD, random_seed=None):
    run_results = {}

    ciw.seed(random_seed)

    model = get_model(experiment)

    sim_engine = ciw.Simulation(model)

    sim_engine.simulate_until_max_time(rc_period)

    recs = sim_engine.get_all_records()

    op_servicetimes = [r.service_time for r in recs if r.node == 1]
    nurse_servicetimes = [r.service_time for r in recs if r.node == 2]

    op_waits = [r.waiting_time for r in recs if r.node == 1]
    nurse_waits = [r.waiting_time for r in recs if r.node == 2]

    run_results["01_mean_waiting_time"] = np.mean(op_waits)
    run_results["02_operator_util"] = (
        sum(op_servicetimes) / (rc_period * experiment.n_operators)
    ) * 100.0
    run_results["03_mean_nurse_waiting_time"] = np.mean(nurse_waits)
    run_results["04_nurse_util"] = (
        sum(nurse_servicetimes) / (rc_period * experiment.n_nurses)
    ) * 100.0

    return run_results, recs


def multiple_replications(experiment, rc_period=RESULTS_COLLECTION_PERIOD, n_reps=5):
    results = []
    logs = []

    for rep in range(n_reps):
        run_result, log = single_run(experiment, rc_period)
        results.append(run_result)
        logs.append(log)

    df_results = pd.DataFrame(results)
    df_results.index = np.arange(1, len(df_results) + 1)
    df_results.index.name = "rep"

    return df_results, logs
N_OPERATORS = 18
N_NURSES = 9
RESULTS_COLLECTION_PERIOD = 1000

user_experiment = Experiment(
    n_operators=N_OPERATORS, n_nurses=N_NURSES, chance_callback=0.4
)

# run multiple replications
results, logs = multiple_replications(user_experiment, n_reps=10)

Here, the ‘logs’ object is the result of running sim_engine.get_all_records()

However, note that while we run multiple replications, we only pass the records for a single replication to the event_log_from_ciw_recs function.

While we’ve done multiple replications, for the purpose of the animation we want only a single set of logs, so we will extract those from the logs variable we created.

# the 'logs' object contains a list, where each entry is the recs object for that run
logs_run_1 = logs[0]

print(len(logs_run_1))
2227

Let’s look at the first row of our result. What do we have?

logs[0][0]
Record(id_number=1533, customer_class='Customer', original_customer_class='Customer', node=1, arrival_date=928.2665472465916, waiting_time=0.4604764322907613, service_start_date=928.7270236788823, service_time=5.280845677909156, service_end_date=934.0078693567915, time_blocked=0.0, exit_date=934.0078693567915, destination=2, queue_size_at_arrival=18, queue_size_at_departure=17, server_id=3, record_type='service')

Let’s print all of the outputs for a single individual.

[print(log) for log in logs_run_1 if log.id_number == 500]
Record(id_number=500, customer_class='Customer', original_customer_class='Customer', node=1, arrival_date=282.0604797501058, waiting_time=0.0, service_start_date=282.0604797501058, service_time=6.668194741388618, service_end_date=288.7286744914944, time_blocked=0.0, exit_date=288.7286744914944, destination=2, queue_size_at_arrival=12, queue_size_at_departure=18, server_id=13, record_type='service')
Record(id_number=500, customer_class='Customer', original_customer_class='Customer', node=2, arrival_date=288.7286744914944, waiting_time=30.463225055971407, service_start_date=319.1918995474658, service_time=18.31678922087241, service_end_date=337.50868876833823, time_blocked=0.0, exit_date=337.50868876833823, destination=-1, queue_size_at_arrival=25, queue_size_at_departure=28, server_id=9, record_type='service')
[None, None]

It looks like we get one entry per node that is visited.

Let’s now use some functions vidigi provides to turn it into a format that can be used for our animation.

The _from_ciw_recs family of functions

Tip

Prior to vidigi 2.0.0, we would have made use of the event_log_from_ciw_recs helper function from vidigi.utils to automatically reshape ciw logs into the correct format for vidigi to work with. This would return the ciw logs as a dataframe in the format that vidigi can use. However, from v2.0.0, the recommendation has changed to using event_logger_from_ciw_recs returning a vidigi EventLogger object, which opens up a wide range of additional visualisations and works more comfortably with the vidigi animation functions. v2.0.0 also introduced the new trial_log_from_ciw_recs, which we will explore as well.

First, we’ll have a go with event_logger_from_ciw_recs. We earlier imported this from the ciw submodule of vidigi with from vidigi.ciw import event_logger_from_ciw_recs.

For each node, we must pass in an appropriate name. Vidigi will use these and append ’_begins’ and ’_ends’, as well as calculating arrivals and departures from the model and creating resource IDs to allow it to correctly show the utilisation of a resource at each step.

# let's now try turning this into an event log
event_log_test = event_logger_from_ciw_recs(
    logs_run_1, node_name_list=["operator", "nurse"]
)

event_log_test
<vidigi.logging.EventLogger at 0x7f75b6f0e5a0>

Let’s use the .to_dataframe() method of our event logger to take a peek at the created dataframe.

event_log_test.to_dataframe().head(25)
entity_id event_type event time pathway resource_id
0 1 arrival_departure arrival 1.169860 Model NaN
1 1 queue operator_wait_begins 1.169860 Model NaN
2 1 resource_use operator_begins 1.169860 Model 1.0
3 1 resource_use_end operator_ends 7.470014 Model 1.0
4 1 arrival_departure depart 7.470014 Model NaN
5 2 arrival_departure arrival 1.709839 Model NaN
6 2 queue operator_wait_begins 1.709839 Model NaN
7 2 resource_use operator_begins 1.709839 Model 2.0
8 2 resource_use_end operator_ends 10.453672 Model 2.0
9 2 arrival_departure depart 10.453672 Model NaN
10 3 arrival_departure arrival 2.103499 Model NaN
11 3 queue operator_wait_begins 2.103499 Model NaN
12 3 resource_use operator_begins 2.103499 Model 3.0
13 3 resource_use_end operator_ends 9.301549 Model 3.0
14 3 queue nurse_wait_begins 9.301549 Model NaN
15 3 resource_use nurse_begins 9.301549 Model 1.0
16 3 resource_use_end nurse_ends 21.038331 Model 1.0
17 3 arrival_departure depart 21.038331 Model NaN
18 4 arrival_departure arrival 3.516870 Model NaN
19 4 queue operator_wait_begins 3.516870 Model NaN
20 4 resource_use operator_begins 3.516870 Model 4.0
21 4 resource_use_end operator_ends 11.508091 Model 4.0
22 4 queue nurse_wait_begins 11.508091 Model NaN
23 4 resource_use nurse_begins 11.508091 Model 2.0
24 4 resource_use_end nurse_ends 30.444661 Model 2.0

Let’s look at using TrialLogger.

trial_log = trial_logger_from_ciw_recs(
    logs, node_name_list=["operator", "nurse"]
)
trial_log
<vidigi.logging.TrialLogger at 0x7f75b6993230>
trial_log.get_log_by_run(run=5, as_df=True).head(25)
entity_id event_type event time pathway run_number resource_id
0 1 arrival_departure arrival 0.447076 Model 5 NaN
1 1 queue operator_wait_begins 0.447076 Model 5 NaN
2 1 resource_use operator_begins 0.447076 Model 5 1.0
3 1 resource_use_end operator_ends 6.542522 Model 5 1.0
4 1 arrival_departure depart 6.542522 Model 5 NaN
5 2 arrival_departure arrival 0.648398 Model 5 NaN
6 2 queue operator_wait_begins 0.648398 Model 5 NaN
7 2 resource_use operator_begins 0.648398 Model 5 2.0
8 2 resource_use_end operator_ends 7.637025 Model 5 2.0
9 2 arrival_departure depart 7.637025 Model 5 NaN
10 3 arrival_departure arrival 1.661917 Model 5 NaN
11 3 queue operator_wait_begins 1.661917 Model 5 NaN
12 3 resource_use operator_begins 1.661917 Model 5 3.0
13 3 resource_use_end operator_ends 10.216653 Model 5 3.0
14 3 arrival_departure depart 10.216653 Model 5 NaN
15 4 arrival_departure arrival 2.011635 Model 5 NaN
16 4 queue operator_wait_begins 2.011635 Model 5 NaN
17 4 resource_use operator_begins 2.011635 Model 5 4.0
18 4 resource_use_end operator_ends 9.416627 Model 5 4.0
19 4 arrival_departure depart 9.416627 Model 5 NaN
20 5 arrival_departure arrival 2.371072 Model 5 NaN
21 5 queue operator_wait_begins 2.371072 Model 5 NaN
22 5 resource_use operator_begins 2.371072 Model 5 5.0
23 5 resource_use_end operator_ends 9.474446 Model 5 5.0
24 5 queue nurse_wait_begins 9.474446 Model 5 NaN

Using the TrialLogger class gets us access to a wide range of convenience plots and metrics in addition to the animation we will create shortly.

trial_log.plot_queue_size(
    event_list=["operator_wait_begins"],
    limit_duration=RESULTS_COLLECTION_PERIOD,
    every_x_time_units=10
    )

Now we need to create a suitable class to pass in the resource numbers to the animation function.

Like with SimPy, we need to tell vidigi where to put each step on our plot. We will refer to the names we used - so as we named our nodes ‘operator’ and ‘nurse’, we will want

  • arrival
  • operator_wait_begins (to show queueing for the operator)
  • operator_begins (to show resource use of the operator)
  • nurse_wait_begins (to show queuing for the nurse after finishing being seen by the operator)
  • nurse_begins (to show resource use of the nurse)

For the _begins steps, which relat to resource use, we will also pass in a name that relates to the number of resources we need, which we defined in our model_params class above.

So, for the operator_begins step, for example, we tell ut to look for n_operators, whch is one of the parameters in our model_params class. We pass the params class into the animation function.

# Create required event_position_df for vidigi animation
event_position_df = create_event_position_df(
    [
        EventPosition(
            event="operator_wait_begins", x=205, y=270, label="Waiting for Operator"
        ),
        EventPosition(
            event="operator_begins",
            x=210,
            y=210,
            resource="n_operators",
            label="Speaking to Operator",
        ),
        EventPosition(
            event="nurse_wait_begins", x=205, y=110, label="Waiting for Nurse"
        ),
        EventPosition(
            event="nurse_begins",
            x=210,
            y=50,
            resource="n_nurses",
            label="Speaking to Nurse",
        ),
        EventPosition(event="depart", x=270, y=10, label="Exit"),
    ]
)

event_position_df
event x y label resource direction flip_icons resource_icon
0 operator_wait_begins 205 270 Waiting for Operator NaN None None None
1 operator_begins 210 210 Speaking to Operator n_operators None None None
2 nurse_wait_begins 205 110 Waiting for Nurse NaN None None None
3 nurse_begins 210 50 Speaking to Nurse n_nurses None None None
4 depart 270 10 Exit NaN None None None

Finally, we can create the animation.

We can access the main animation function animate_activity_log directly from our TrialLogger object (or an EventLogger object if we’d just done a single run).

We will pass in our resource counts as a dict to the scenario parameter.

trial_log.animate_activity_log(
    run_number=5,
    event_position_df=event_position_df,
    scenario={'n_operators':N_OPERATORS, 'n_nurses': N_NURSES},
    debug_mode=True,
    setup_mode=False,
    every_x_time_units=1,
    include_play_button=True,
    entity_icon_size=20,
    gap_between_entities=8,
    gap_between_queue_rows=25,
    plotly_height=700,
    frame_duration=200,
    plotly_width=1200,
    override_x_max=300,
    override_y_max=300,
    limit_duration=RESULTS_COLLECTION_PERIOD,
    wrap_queues_at=25,
    wrap_resources_at=50,
    step_snapshot_max=75,
    time_display_units="dhm",
    display_stage_labels=True,
)
Animation function called at 13:54:23
Iteration through time-unit-by-time-unit logs complete 13:54:31
Snapshot df concatenation complete at 13:54:31
Reshaped animation dataframe finished construction at 13:54:32
Placement dataframe started construction at 13:54:32
Placement dataframe finished construction at 13:54:32
Output animation generation complete at 13:54:37
Total Time Elapsed: 14.16 seconds
Back to top