Chapter 2: Simulation Runs and Replications

How runSim and parallelProcessJoblib run seeded replications of the air ambulance model and collate their output files.

Written for: Intermediate Coder Confident Coder Decision-makers Domain experts General public

In Streamlit App and Scenario Parameters we saw the app gather a scenario (number of runs, duration, start date, demand multiplier, master seed) and then call into the package to run it. This chapter looks at what is on the other end of that call: the functions in air_ambulance_des/des_parallel_process.py that turn those numbers into many independent runs of the model and a single results file.

Why one run is not enough

Our central use case is the one the app’s “Run” button performs. A user asks for 12 replications of a two-year simulation, starting on 1 January 2023, with the master seed 42. They want a combined event log they can summarise (see Trial Results and Historic Comparison), and they want to be able to run the same scenario tomorrow and get exactly the same numbers.

Those two wants pull in different directions. A discrete event simulation is stochastic: call numbers, patient types and activity durations are all sampled, so one run is just one possible two years. You need several replications to see the spread. Yet the answer must also be reproducible, and the replications must not be copies of each other. The code solves this by deriving one independent seed per replication from a single master seed, then running the replications across CPU cores.

Replication: one complete run of the simulation with its own random seed. Results are compared across replications to see how much chance matters.

Master seed: the single integer the user controls. Every per-run seed is derived from it, so the same master seed always gives the same set of runs.

The pieces

Four functions do the work, plus the model’s own run() method:

Piece	Role
parallelProcessJoblib	Derives per-run seeds and fans the runs out across cores
runSim	Builds and runs one DES_HEMS model, returns its results
write_run_params	Records the settings used, once, in run_params_used.csv
collateRunResults	Merges the per-run CSVs into run_results.csv
DES_HEMS.run	Runs the SimPy clock and writes the run’s CSV

There are two time conventions to keep in mind. The model’s clock is in simulated minutes, so a two-year duration is 2 * 365 * 24 * 60. And the warm-up period is an initial stretch of simulated time that is run but not recorded, so that results are not distorted by the model starting from an empty system.

Warm-up period: simulated minutes at the start of a run during which events happen but no rows are written to the event log.

Using it

The command-line equivalent of the app’s run button is run_parallel_process.py in the repository root. It is a good place to see the use case in its plainest form.

run_parallel_process.py:1-19

from air_ambulance_des.des_parallel_process import (
    removeExistingResults,
    parallelProcessJoblib,
)
from datetime import datetime

if __name__ == "__main__":
    removeExistingResults()
    # parallelProcessJoblib(1, (1*365*24*60), (0*60), datetime.strptime("2023-01-01 05:00:00", "%Y-%m-%d %H:%M:%S"), False, False, 1.0, 1.0, True)
    parallelProcessJoblib(
        12,
        (2 * 365 * 24 * 60),
        (0 * 60),
        datetime.strptime("2023-01-01 05:00:00", "%Y-%m-%d %H:%M:%S"),
        False,
        False,
        1.0,
        1.0,
    )

removeExistingResults() clears old output_run_*.csv and all_results.csv files from data/ first. The call then passes, positionally, 12 runs, a duration of two years in minutes, zero warm-up, the start datetime, amb_data=False (do not model ambulance-service cases alongside HEMS; this is the normal setting), save_params_csv=False (do not write run_params_used.csv from this call), and demand and activity-duration multipliers of 1.0. The master seed is not given, so the default of 42 applies.

The app makes the same call, shown in full in chapter 1, with keyword arguments. The lines that matter here are the master_seed= argument, which comes from session state, and the collateRunResults() call that follows parallelProcessJoblib. The other values are converted as set up in that chapter: duration is days multiplied by 24 * 60, warm-up is hours multiplied by 60, and the demand multiplier is a percentage divided by 100.

WarningThe command-line script does not collate

run_parallel_process.py calls removeExistingResults() and parallelProcessJoblib but never collateRunResults(). Run on its own, it leaves the per-run files rather than a combined run_results.csv. The app does call the collation step, as shown above. If you use the script, call collateRunResults() yourself afterwards.

What happens when you call it

1. parallelProcessJoblib builds a SeedSequence from the master seed and spawns one child per replication, then turns each child into an integer seed.
2. joblib starts up to n_cores worker processes (all cores by default) and calls runSim once per replication, handing each its own seed.
3. Each runSim creates a DES_HEMS model with that seed and calls its run() method. The model simulates the full period and writes output_run_N.csv.
4. The first run (run == 0) also writes run_params_used.csv.
5. When all workers have finished, the caller runs collateRunResults(), which concatenates the per-run files into run_results.csv and deletes them.

sequenceDiagram
    participant App as App or script
    participant PP as parallelProcessJoblib
    participant RS as runSim (per worker)
    participant M as DES_HEMS
    participant D as data folder
    App->>PP: total_runs=3, master_seed=42
    PP->>PP: spawn 3 seeds from master seed
    PP->>RS: run=0, 1, 2 with seeds (in parallel)
    RS->>M: DES_HEMS(run_number, random_seed, ...)
    RS->>M: run()
    M->>D: output_run_N.csv
    RS->>D: run_params_used.csv (run 0 only)
    RS-->>PP: results DataFrame
    PP-->>App: list of DataFrames
    App->>D: collateRunResults() merges into run_results.csv

Figure 1: From one call to a combined results file (three replications shown)

Figure 1 shows that two things come back from a run: a DataFrame returned to the caller, and a CSV on disk. The app ignores the returned list and relies on the files.

Under the hood: deriving the seeds

The seeding is the most important design decision here, so it comes first.

air_ambulance_des/des_parallel_process.py:294-319

    # Generate a number of uncorrelated seeds that will always be the same given the same
    # master seed (which is determined as a parameter)
    # We start with a SeedSequence from the master seed, and then generate a number of
    # child SeedSequences equal to the total number of runs
    seed_sequence = SeedSequence(master_seed).spawn(total_runs)
    # We then turn these seeds into integer random numbers, and we will pass a different seed
    # into each run of the simulation.
    seeds = [i.generate_state(1)[0] for i in seed_sequence]

    # Run the simulation in parallel, using all available cores
    return Parallel(n_jobs=n_cores)(
        delayed(runSim)(
            run=run,
            total_runs=total_runs,
            sim_duration=sim_duration,
            warm_up_time=warm_up_time,
            sim_start_date=sim_start_date,
            amb_data=amb_data,
            random_seed=seeds[run],
            save_params_csv=save_params_csv,
            demand_increase_percent=demand_increase_percent,
            activity_duration_multiplier=activity_duration_multiplier,
            print_debug_messages=print_debug_messages,
        )
        for run in range(total_runs)
    )

SeedSequence(master_seed).spawn(total_runs) is NumPy’s mechanism for producing statistically independent child sequences from one parent. Simply using seeds 42, 43, 44, ... for the replications would work mechanically, but neighbouring integer seeds are not guaranteed to give uncorrelated streams. Spawning avoids that. generate_state(1)[0] then collapses each child into a single integer, because runSim and DES_HEMS take an int seed.

Because the seeds are computed up front, in the parent process, replication n gets the same seed whatever order the workers happen to finish in, and whatever the number of cores. That is what makes a run reproducible. The Parallel(n_jobs=n_cores)(delayed(runSim)(...) for run in range(total_runs)) pattern is joblib’s standard way of mapping a function over a range of arguments in separate processes.

The seeding does not stop there. Inside the model, the per-run seed is spawned again, this time into many named random streams, one per sampling purpose. That is the subject of Sampling Distributions and Random Streams. The start of DES_HEMS.__init__ is where the hand-over happens:

air_ambulance_des/des_hems.py:56-71

        self.run_number = run_number + 1  # Add 1 so we don't have a run_number 0
        self.sim_duration = sim_duration
        self.warm_up_duration = warm_up_duration
        self.sim_start_date = sim_start_date

        self.random_seed = random_seed
        self.random_seed_sequence = SeedSequence(self.random_seed)

        self.print_debug_messages = print_debug_messages

        self.utils = Utils(
            master_seed=self.random_seed_sequence.spawn(1)[0],
            print_debug_messages=self.print_debug_messages,
        )

        self.utils.setup_seeds()

Note self.run_number = run_number + 1. The run argument counts from 0 (which is how runSim recognises “the first run”), but the run number used in file names starts from 1, so the first replication writes output_run_1.csv.

Under the hood: one run

runSim is a thin wrapper. Most of its length is docstring and progress messages; the work is constructing the model, running it and, for the first replication, saving the parameters.

air_ambulance_des/des_parallel_process.py:167-197

    start = time.process_time()

    print(f"{Utils.current_time()}: Run {run + 1} of {total_runs}")
    logging.debug(f"{Utils.current_time()}: Run {run + 1} of {total_runs}")

    # print(f"Sim start date is {sim_start_date}")
    daa_model = DES_HEMS(
        run_number=run,
        sim_duration=sim_duration,
        warm_up_duration=warm_up_time,
        sim_start_date=sim_start_date,
        amb_data=amb_data,
        random_seed=random_seed,
        demand_increase_percent=demand_increase_percent,
        activity_duration_multiplier=activity_duration_multiplier,
        print_debug_messages=print_debug_messages,
    )
    daa_model.run()

    print(
        f"{Utils.current_time()}: Run {run + 1} took {round((time.process_time() - start) / 60, 1)} minutes to run"
    )
    logging.debug(
        f"{Utils.current_time()}: Run {run + 1} took {round((time.process_time() - start) / 60, 1)} minutes to run"
    )

    # SR NOTE: This could cause issues if we decide to use 1 as the starting number of runs
    if (run == 0) and (save_params_csv):
        write_run_params(daa_model)

    return daa_model.results_df

Three details are worth noticing.

• time.process_time() measures CPU time used by the current process. Because each replication runs in its own worker process, this is a fair per-run timing, though it is not wall-clock time.
• The run == 0 condition is why the parameters file is written exactly once. All replications share the same settings, so one copy is enough. The SR NOTE comment in the source points out that this would break if run numbering ever started at 1.
• write_run_params(daa_model) reads the settings back from the model (duration, warm-up, start date, multipliers) and writes them to run_params_used.csv, including derived values such as the end date and the warm-up end date. Trial Results and Historic Comparison reads this file back later.

air_ambulance_des/des_parallel_process.py:62-97

    # Ensure sim_start_date is a datetime object
    sim_start_date = model.sim_start_date
    if isinstance(sim_start_date, str):
        sim_start_date = datetime.fromisoformat(
            sim_start_date
        )  # Convert string to datetime

    sim_end_date = sim_start_date + timedelta(minutes=model.sim_duration)
    warm_up_end_date = sim_start_date + timedelta(minutes=model.warm_up_duration)

    params_df = pd.DataFrame.from_dict(
        {
            "sim_duration": [model.sim_duration],  # model time unit is minutes
            "sim_duration_hours": [model.sim_duration / 60],
            "sim_duration_days": [model.sim_duration / 60 / 24],
            "warm_up_duration": [model.warm_up_duration],
            "sim_start_date": [sim_start_date],
            "sim_end_date": [sim_end_date],
            "warm_up_end_date": [warm_up_end_date],
            "amb_data": [model.amb_data],
            "model_exec_time": [datetime.now()],
            # Assuming summer hours are quarters 2 and 3 i.e. April-September
            # This is defined in class_hems and will need updating here too
            "summer_start_date": [f"{sim_start_date.year}-04-01"],
            "winter_start_date": [f"{sim_start_date.year}-10-01"],
            "activity_duration_multiplier": [model.activity_duration_multiplier],
        },
        orient="index",
        columns=["value"],
    )

    params_df.index.name = "parameter"

    params_df.to_csv(
        f"{Utils.RESULTS_FOLDER}/run_params_used.csv", header="column_names"
    )

The summer and winter start dates are hard-coded as 1 April and 1 October of the start year. The comment in the code notes that this duplicates a definition in class_hems and must be kept in step with it.

Under the hood: inside DES_HEMS.run

The model’s run() method is short, because the real work lives in processes registered on the SimPy environment (covered in Calls and Patient Journeys).

air_ambulance_des/des_hems.py:812-832

    def run(self) -> None:
        """
        Function to start the simulation.

        """
        self.debug(
            f"HEMS class initialised with the following: run {self.run_number}, duration {self.sim_duration}, warm-up {self.warm_up_duration}, start date {self.sim_start_date}, demand increase multiplier {self.demand_increase_percent}, activity duration multiplier {self.activity_duration_multiplier}"
        )

        # Start entity generators
        self.env.process(self.generate_calls())

        # Run simulation
        self.env.run(until=(self.sim_duration + self.warm_up_duration))

        # Convert results to a dataframe
        self.convert_results_to_df(results=self.results_list)

        # Write run results to file
        self.write_all_results()
        self.write_run_results()

The simulation runs until sim_duration + warm_up_duration: the warm-up is added on top of the requested duration, not carved out of it. A one-year request with 24 hours of warm-up simulates 365 days plus one day, and records the last 365.

The exclusion of the warm-up happens when events are logged, not here. Every call to add_patient_result_row checks the clock before storing a row:

air_ambulance_des/des_hems.py:764-765

        if self.env.now >= self.warm_up_duration:
            self.store_patient_results(results)

So during the warm-up, patients still arrive and still occupy helicopters and cars, but nothing is added to results_list. By the time recording begins, the model is already in a realistic state.

Finally, run() writes two files. They are worth distinguishing:

air_ambulance_des/des_hems.py:782-810

    def write_all_results(self) -> None:
        """
        Writes the content of `result_df` to a csv file
        """
        # https://stackoverflow.com/a/30991707/3650230

        # Check if file exists...if it does, append data, otherwise create a new file with headers matching
        # the column names of `results_df`
        if not os.path.isfile(self.all_results_location):
            self.results_df.to_csv(self.all_results_location, header="column_names")
        else:  # else it exists so append without writing the header
            self.results_df.to_csv(self.all_results_location, mode="a", header=False)

    def write_run_results(self) -> None:
        """
        Writes the content of `result_dfs` to csv files that contains only the results from the
        single run

        Note that this cannot be done in a similar manner to write_all_results due to the impacts of
        the parallisation approach that is taken with joblib - depending on process timings it can lead to
        not all
        """

        self.results_df.to_csv(
            f"{Utils.RESULTS_FOLDER}/output_run_{self.run_number}.csv",
            header="column_names",
            index=True,
            encoding="utf-8-sig",
        )

write_all_results appends to a single shared all_results.csv, which is why removeExistingResults deletes it before a fresh batch. write_run_results writes output_run_N.csv, one file per replication. The docstring explains why the combined file is not simply the shared one: with joblib, several workers could be appending at once, and depending on timing not all rows would be reliably written. Separate per-run files cannot collide, so they are the safe route to a combined result.

Under the hood: collation

air_ambulance_des/des_parallel_process.py:200-211

def collateRunResults() -> None:
    """
    Collates results from a series of runs into a single csv
    """
    matching_files = glob.glob(os.path.join(Utils.RESULTS_FOLDER, "output_run_*.csv"))

    combined_df = pd.concat([pd.read_csv(f) for f in matching_files], ignore_index=True)

    combined_df.to_csv(Utils.RUN_RESULTS_CSV, index=True, encoding="utf-8-sig")

    for file in matching_files:
        os.remove(file)

collateRunResults finds every output_run_*.csv in Utils.RESULTS_FOLDER (the data/ folder), reads them with pandas, concatenates them with a fresh index and writes Utils.RUN_RESULTS_CSV. It then deletes the per-run files. Each row already carries its run number, so the combined file still distinguishes replications.

This is the run_results.csv that the results classes consume. Its sibling, run_params_used.csv, describes the settings under which it was produced.

NoteRefresher: why joblib uses processes

Python threads share one interpreter lock, so CPU-bound work such as a simulation gains little from them. joblib’s default backend starts separate worker processes, each with its own interpreter and memory. That is why the functions return DataFrames or write files rather than updating shared variables, and why the seeds are computed in the parent and passed in as plain integers.

Conclusion

A scenario becomes a batch of replications through three steps: parallelProcessJoblib spawns one independent seed per replication from the master seed, runSim builds and runs a DES_HEMS model for each, and collateRunResults merges the per-run CSVs into run_results.csv. Warm-up time is simulated but not logged, and run_params_used.csv records the settings alongside.

Each model is handed an integer seed, but we have not yet seen what it does with it. Next, in Sampling Distributions and Random Streams, we’ll see how that seed is split into named streams and how fitted distributions are sampled. For the project background, see the project overview.