flowchart TD
M["Master seed (chapter 2)"] --> R["Per-run seed"]
R --> D["DES_HEMS: SeedSequence of the run seed"]
D --> U["spawn(1): master_seed for Utils"]
U --> S["setup_seeds: spawn one child per module key"]
S --> K1["rngs calls_per_day"]
S --> K2["rngs ampds_code_selection"]
S --> K3["rngs ... 22 named streams"]
U --> B["build_seeded_distributions: spawn one child per entry in activity_time_distributions.txt"]
B --> A1["SeededDistribution (car, time_on_scene)"]
B --> A2["SeededDistribution (helicopter, time_to_scene)"]
Chapter 3: Sampling Distributions and Random Streams
How Utils loads the fitted distributions and samples from them with named, independently seeded random streams.
In Simulation Runs and Replications we saw that a master seed is split into one independent seed per replication, and that each replication builds its own DES_HEMS model from its seed. This chapter looks at what the model does with that seed. Every random draw in a run (how many calls arrive today, which AMPDS card a call gets, how long a helicopter takes to reach the scene) comes from a named random stream owned by one class, Utils, in air_ambulance_des/utils.py. An AMPDS card is the triage code the dispatch system assigns to a call; it is defined in Calls and Patient Journeys.
The problem: random, but repeatable
Our central use case is this. A user runs the same scenario twice with the same master seed and expects identical results. Then they change one thing, say the activity-time distribution for cars, and expect the only differences to come from that change. They do not want a different number of calls per day appearing as a side effect.
A single shared random number generator cannot give you that. If the car’s time-on-scene sampler draws one extra number, every later draw from the shared generator shifts, and the whole run changes. Utils avoids this by giving each sampling purpose its own named random stream.
Think of one shared deck of cards versus a separate sealed deck for each purpose. If you draw an extra card for car timings from the shared deck, every later draw in the run changes. With a deck per purpose, only that purpose’s sequence moves.
Named random stream: a separate numpy random generator reserved for one sampling purpose, such as rngs["calls_per_day"]. See the glossary.
The chapter follows that use case in three steps: where Utils gets its distributions, how it creates the streams, and how individual selection methods use them.
What Utils loads
Utils is a plain class with no SimPy dependency. Its class attributes hold the folder locations; everything else is read in __init__.
air_ambulance_des/utils.py:34-45
PACKAGE_ROOT = Path(__file__).resolve().parent # points to air_ambulance_des/
REPO_ROOT = PACKAGE_ROOT.parent # points to repo root
class Utils:
RESULTS_FOLDER = REPO_ROOT / "data"
ACTUAL_DATA_FOLDER = REPO_ROOT / "actual_data"
HISTORICAL_FOLDER = REPO_ROOT / "historical_data"
DISTRIBUTION_FOLDER = REPO_ROOT / "distribution_data"
ALL_RESULTS_CSV = RESULTS_FOLDER / "all_results.csv"
RUN_RESULTS_CSV = RESULTS_FOLDER / "run_results.csv"PACKAGE_ROOT and REPO_ROOT are worked out from the file’s own location, so the paths do not depend on the directory you launch Python from. Three folders matter here:
DISTRIBUTION_FOLDER(distribution_data/) holds the fitted inputs, produced offline by the pipeline in Distribution Fitting. Nothing in a simulation run writes to it.ACTUAL_DATA_FOLDER(actual_data/) holds the user-editable inputs from Streamlit App and Scenario Parameters.RESULTS_FOLDER(data/) is where the run outputs from chapter 2 land.
distribution_data folder: the repository folder holding the CSV, text and JSON files of fitted probabilities and distributions that the simulation samples from.
The constructor then reads those files into attributes. They fall into three shapes.
Probability tables are CSVs read into DataFrames, for example hour_by_ampds_card_probs.csv (probability of each AMPDS card by hour), enhanced_or_critical_care_by_ampds_card_probs.csv, vehicle_type_by_quarter_probs.csv and the patient outcome and HEMS result tables. Each has a proportion column that sums to one within a group, which is what makes weighted sampling possible later.
Fitted distributions are stored as Python literals in text files and parsed with ast.literal_eval:
air_ambulance_des/utils.py:147-154
# Read in activity time distribution data into a dictionary
activity_time_data = []
with open(
self.DISTRIBUTION_FOLDER / "activity_time_distributions.txt", "r"
) as inFile:
activity_time_data = ast.literal_eval(inFile.read())
inFile.close()
self.activity_time_distr = activity_time_dataEach entry in activity_time_distributions.txt names a HEMS result, a vehicle type, a time type and a best-fit scipy distribution with its parameters, for example a Cauchy with loc and scale for helicopter time_allocation. At this point self.activity_time_distr is still a list of raw dictionaries. It is turned into something sampleable in build_seeded_distributions(), below.
Daily call counts are not a fitted curve at all. inc_per_day_samples.json holds lists of observed daily call counts, keyed by season (winter, summer) and by quarter (Q1 to Q4). Sampling a day’s demand means picking one of those historical days at random.
Finally, Utils loads the upper allowable time bounds and reshapes them for fast lookup:
air_ambulance_des/utils.py:125-128
# Import maximum call duration times
self.min_max_values_df = pd.read_csv(
self.ACTUAL_DATA_FOLDER / "upper_allowable_time_bounds.csv"
)air_ambulance_des/utils.py:177-182
# Turn the min-max activity times data into a format that supports easier/faster
# lookups
self.min_max_cache = {
row["time"]: (row["min_value_mins"], row["max_value_mins"])
for _, row in self.min_max_values_df.iterrows()
}Upper allowable time bounds: per-activity minimum and maximum durations in minutes, read from actual_data/upper_allowable_time_bounds.csv. Sampled activity times outside them are rejected.
min_max_cache is a dictionary from a time type such as time_on_scene to a (min, max) tuple. The CSV rows look like time_on_scene,1.0,119.0. Because these bounds live in actual_data/, a user can change them as part of a scenario, unlike the fitted files.
One more file is loaded defensively: the table for ad hoc unavailability (ad_hoc_unavailability.csv) sits inside a try/except FileNotFoundError, because the fitting pipeline itself creates Utils before that file exists.
Creating the named streams
Loading does not touch randomness. Seeding happens in a separate call, setup_seeds(). First, here is how DES_HEMS builds and seeds its Utils:
air_ambulance_des/des_hems.py:61-71
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()The random_seed is the per-run seed that runSim received, so everything below is a deterministic function of it. DES_HEMS wraps it in a SeedSequence, spawns one child to hand to Utils as its master_seed, and calls setup_seeds().
SeedSequence spawning: numpy’s way of deriving statistically independent child seeds from a parent. Each spawn() call returns new children and remembers how many it has already produced.
air_ambulance_des/utils.py:184-221
def setup_seeds(self):
#########################################################
# Control generation of required random number streams #
#########################################################
# Spawn substreams for major simulation modules
module_keys = [
"activity_times",
"ampds_code_selection",
"care_category_selection",
"callsign_group_selection",
"vehicle_type_selection",
"hems_result_by_callsign_group_and_vehicle_type_selection",
"hems_result_by_care_category_and_helicopter_benefit_selection",
"hems_results_by_patient_outcome_and_quarter_and_vehicle_type_and_callsign_group_selection",
"pt_outcome_selection",
"sex_selection",
"age_sampling",
"calls_per_day",
"calls_per_hour",
"predetermine_call_arrival",
"call_iat",
"helicopter_benefit_from_reg",
"hems_case",
"ad_hoc_reason_selection",
"know_heli_benefit",
"helicopter_benefit_from_cc",
"helicopter_benefit_from_ec",
"know_cc_ec_benefit",
]
# Efficiently spawn substreams
spawned = self.master_seed_sequence.spawn(len(module_keys))
# Map substreams and RNGs to keys
self.seed_substreams = dict(zip(module_keys, spawned))
self.rngs = {key: default_rng(ss) for key, ss in self.seed_substreams.items()}
self.build_seeded_distributions(self.master_seed_sequence)The method works in four moves:
module_keysis a hard-coded list of every sampling purpose, one name each. Some are obvious (ampds_code_selection,calls_per_day), some are more specific (helicopter_benefit_from_cc,know_heli_benefit).self.master_seed_sequence.spawn(len(module_keys))creates one childSeedSequenceper name in a single call.dict(zip(...))pairs each name with its child, anddefault_rng(ss)turns each child into anumpy.random.Generator. The result isself.rngs, a dictionary of independent generators.build_seeded_distributions()is then called with the same sequence, so the activity-time distributions get their own streams too (see below).
Why this arrangement works: children of one SeedSequence are independent of each other, and each depends only on the parent seed and its position in the list. Drawing a thousand numbers from rngs["ampds_code_selection"] does not change what rngs["calls_per_day"] produces next. That is exactly the property the use case needs. (Two of the named streams, activity_times and predetermine_call_arrival, are reserved but never drawn from: activity times use the per-distribution generators described below.)
Figure 1 shows the whole seed tree from master seed down to individual streams.
WarningEditing the key list shifts later streams
A child’s stream depends on its position in the spawn order, and spawn() is stateful. Appending a name to the end of module_keys leaves the earlier named streams unchanged, but inserting one in the middle moves every stream after it. Adding keys at the end also advances the spawn counter, so the activity-time generators created afterwards by build_seeded_distributions() receive different seeds. Tests such as tests/test_reproducibility.py are the place to check that a change like this does what you intended.
Sampling a single value
With streams in place, the selection methods are short. The most common shape is a weighted choice from a probability table.
air_ambulance_des/utils.py:289-301
def ampds_code_selection(self, hour: int) -> int:
"""
This function will allocate and return an AMPDS card category
based on the hour of day
"""
df = self.hour_by_ampds_df[self.hour_by_ampds_df["hour"] == hour]
return pd.Series.sample(
df["ampds_card"],
weights=df["proportion"],
random_state=self.rngs["ampds_code_selection"],
).iloc[0]ampds_code_selection() filters the table to the current hour, then uses pandas.Series.sample with the proportion column as weights. The important detail is random_state=self.rngs["ampds_code_selection"]: pandas accepts a Generator, so the draw consumes that named stream and no other. care_category_selection(), callsign_group_selection(), vehicle_type_selection_qtr() and pt_outcome_selection() follow the same pattern with their own stream names. How the patient generator chains these calls together is the subject of Calls and Patient Journeys.
NoteRefresher: weighted sampling with pandas
Series.sample(weights=w) normalises w and draws rows with those probabilities. With random_state set to a numpy Generator rather than an integer, repeated calls continue that generator’s sequence instead of restarting it. Passing an integer would give the same answer every time.
Daily call counts
Demand for the day is sampled by inc_per_day_samples():
air_ambulance_des/utils.py:652-671
def inc_per_day_samples(
self, quarter: int, quarter_or_season: str = "season"
) -> float:
"""
This function will return a dictionary containing
the distribution and parameters for the distribution
that match the provided HEMS vehicle type and time type
"""
season = "summer" if quarter in [2, 3] else "winter"
if quarter_or_season == "season":
return self.rngs["calls_per_day"].choice(
self.incident_per_day_samples[season]
)
else:
return self.rngs["calls_per_day"].choice(
self.incident_per_day_samples[f"Q{quarter}"]
)The quarter_or_season argument chooses between the winter/summer lists and the Q1..Q4 lists. DES_HEMS sets it to "quarter" (daily_calls_by_quarter_or_season in des_hems.py), then multiplies the chosen count by the demand multiplier and truncates it to an integer. There is also an older inc_per_day() that samples a fitted distribution instead; the model’s call site for it is commented out in favour of the empirical samples.
Activity times
Durations are where the bounds and the per-distribution streams come together.
air_ambulance_des/utils.py:580-607
def activity_time(self, vehicle_type: str, time_type: str) -> float:
"""
This function will return a dictionary containing
the distribution and parameters for the distribution
that match the provided HEMS vehicle type and time type
"""
dist = self.activity_time_distr.get((vehicle_type, time_type))
# print(dist)
if dist is None:
raise ValueError(f"No distribution found for ({vehicle_type}, {time_type})")
try:
min_time, max_time = self.min_max_cache[time_type]
# print(f"Min time {time_type}: {min_time}")
# print(f"Max time {time_type}: {max_time}")
except KeyError:
raise ValueError(f"Min/max bounds not found for time_type='{time_type}'")
sampled_time = -1
while not (min_time <= sampled_time <= max_time):
sampled_time = dist.sample()
# print(sampled_time)
return sampled_timeActivity time: the sampled duration, in simulated minutes, of one stage of a job (for example time_on_scene) for a given vehicle type.
activity_time(vehicle_type, time_type) looks up the distribution for that pair, looks up the bounds for the time type, and then rejection samples: it keeps drawing until the value lies inside [min_time, max_time]. This is why a rare tail value from a Cauchy fit cannot produce a two-day on-scene time. A missing pair or missing bounds raises a ValueError rather than silently returning a default.
dist here is not a raw dictionary. It is a SeededDistribution, which build_seeded_distributions() put in self.activity_time_distr in place of the loaded list.
SeededDistribution
air_ambulance_des/utils.py:1053-1059
class SeededDistribution:
def __init__(self, dist, rng, **kwargs):
self.dist = dist(**kwargs)
self.rng = rng
def sample(self, size=None):
return self.dist.rvs(size=size, random_state=self.rng)SeededDistribution: a small wrapper pairing a frozen scipy distribution with its own numpy generator, so sample() needs no random-state argument.
A frozen scipy distribution (dist(**kwargs)) fixes its parameters. SeededDistribution adds the generator, and sample() passes it as random_state to rvs. The result is an object you can call .sample() on and trust to be reproducible, which is what activity_time() does. build_seeded_distributions() creates one per fitted entry:
air_ambulance_des/utils.py:902-903
n = len(self.activity_time_distr)
rngs = [default_rng(s) for s in seed_seq.spawn(n)]air_ambulance_des/utils.py:928-955
seeded_distributions = {}
# print(activity_time_distr)
# print(len(self.activity_time_distr))
i = 0
for entry, rng in zip(self.activity_time_distr, rngs):
i += 1
vt = entry["vehicle_type"]
tt = entry["time_type"]
best_fit = entry["best_fit"]
# dist_name = best_fit['dist'].lower()
dist_name = list(best_fit.keys())[0]
dist_cls = dist_map.get(dist_name)
self.debug(
f"{i}/{len(self.activity_time_distr)} for {vt} {tt} set up {dist_name} {dist_cls} with params: {best_fit[dist_name]}"
)
if dist_cls is None:
raise ValueError(f"Unsupported distribution type: {dist_name}")
params = best_fit[dist_name]
seeded_distributions[(vt, tt)] = SeededDistribution(dist_cls, rng, **params)
self.activity_time_distr = seeded_distributions
return seeded_distributionsThe first block spawns one generator per fitted entry. Because the fitting pipeline writes one entry per HEMS result as well as per vehicle type and time type, the file holds more entries than there are (vehicle type, time type) pairs (50 entries for 12 pairs in the shipped file). The dictionary assignment keeps only the last entry for each pair, so the HEMS-result stratification is effectively ignored at sampling time and earlier generators are discarded (see Distribution Fitting). The loop then reads the distribution’s name from the single key of best_fit, looks the class up in a dist_map (defined just above, mapping names such as "gamma", "lognorm" or "cauchy" to scipy classes), and builds the wrapper. An unsupported name raises ValueError. The final two lines replace self.activity_time_distr with a dictionary keyed by (vehicle_type, time_type), which is the key activity_time() looks up.
Warningsetup_seeds() is not idempotent
build_seeded_distributions() overwrites self.activity_time_distr with the tuple-keyed dictionary. Calling setup_seeds() a second time on the same Utils would then iterate over that dictionary’s keys and fail on entry["vehicle_type"]. In normal use DES_HEMS calls it exactly once per run, so this never bites, but it matters if you write a test that reuses a Utils object.
Fitted-distribution sampling for age
Not every fitted distribution is wrapped. Ages use sample_from_distribution(), which takes a {name: params} dictionary and an explicit generator:
air_ambulance_des/utils.py:689-698
if len(distr) != 1:
raise ValueError("Expected one distribution name in distr dictionary")
dist_name, params = list(distr.items())[0]
sci_distr = getattr(scipy.stats, dist_name)
while True:
sampled_value = np.floor(sci_distr.rvs(random_state=rng, **params))
if sampled_value > 0:
return sampled_valueIt resolves the scipy distribution by name with getattr(scipy.stats, dist_name), draws with random_state=rng, floors the value, and loops until it is positive. age_sampling() calls it with rngs["age_sampling"] and repeats until the age is no more than max_age.
Note
sex_selection() and age_sampling() carry a # TODO: RANDOM SEED SETTING comment, but both already draw from self.rngs["sex_selection"] and self.rngs["age_sampling"]. The comment looks like a leftover from before the named streams were introduced.
Ad hoc unavailability
Ad hoc unavailability: a resource that is on shift and serviced but still cannot attend, for a reason such as weather or crew. The probabilities come from historical data.
sample_ad_hoc_reason(hour, quarter, registration) is a different kind of sampler: it returns a reason string such as "available", "weather" or "crew", and it is called when resources are being allocated (see Resource Allocation and HEMSAvailability). The table ad_hoc_unavailability.csv has one row per registration, six-hour bin, quarter and reason. The method first maps the hour onto a bin (00-05, 06-11, 12-17, 18-23) and filters the table. If that registration has no rows, it falls back to every registration (the detail is in the collapsed box below).
NoteDeep dive: how missing ad hoc data is filled in
air_ambulance_des/utils.py:979-992
if subset.empty:
self.debug(
f"[AD HOC] No data for {registration} in bin {bin_label}, Q{quarter}. Falling back to all registrations."
)
subset = self.ad_hoc_probs[
(self.ad_hoc_probs["six_hour_bin"] == bin_label)
& (self.ad_hoc_probs["quarter"] == quarter)
]
if subset.empty:
self.debug(
f"[AD HOC] No data at all for bin {bin_label}, Q{quarter}. Defaulting to 'available'."
)
return "available"
Missing probabilities are then imputed, first from the registration’s average across bins in that quarter, then from the bin’s average across registrations. If anything is still missing, or the total is zero, the method gives up and returns "available". Each fallback logs through self.debug(), which is silent unless print_debug_messages was set.
If the data is still thin after those fallbacks, the method simply returns "available". Once the probabilities are usable the draw is the same pattern as before:
air_ambulance_des/utils.py:1044-1050
norm_probs = filled_subset["probability"] / total_prob
sampled_reason = self.rngs["ad_hoc_reason_selection"].choice(
filled_subset["reason"].tolist(), p=norm_probs.tolist()
)
return sampled_reasonThe probabilities are renormalised so they sum to one, then rngs["ad_hoc_reason_selection"] picks a reason. Defaulting to "available" whenever the data is thin is a deliberate-looking choice: missing data never makes a vehicle spuriously unavailable.
Putting it together
Figure 2 shows a run being set up and then sampling one activity time.
sequenceDiagram
participant RS as runSim
participant DH as DES_HEMS
participant U as Utils
participant SD as SeededDistribution
RS->>DH: create with run seed
DH->>U: Utils(master_seed=child of SeedSequence)
U->>U: __init__ loads distribution_data and bounds
DH->>U: setup_seeds()
U->>U: spawn named rngs, build SeededDistributions
Note over DH,U: simulation runs
DH->>U: activity_time("car", "time_on_scene")
U->>SD: sample()
SD-->>U: value (rejected until within bounds)
U-->>DH: minutes
The flow has a nice property for debugging. Because every stream is derived from the run seed alone, you can construct a DES_HEMS (or just a Utils) with the same seed in a test and call a single selection method to see exactly the value the full run would have produced for that stream’s first draw.
Conclusion
Utils reads the fitted files from distribution_data/ and the time bounds from actual_data/, then setup_seeds() gives each sampling purpose its own generator via SeedSequence.spawn. Selection methods draw weighted choices or empirical samples from those generators, and activity times come from SeededDistribution objects with rejection against the allowable bounds. Together this makes a run repeatable and keeps changes in one part of the model from disturbing the random numbers used elsewhere.
Next, in Calls and Patient Journeys, we’ll see DES_HEMS call these methods: sampling the day’s call count, creating each Patient, and stepping it through the stages of a job.