Benchmark: An Arrayed Model
benchmark, performance, rust, execution engine, arrays, multidimensional, system dynamics, bptk, bptk-py, python
Benchmark: An Arrayed Model
A scalar model grows in one direction: more timesteps. An arrayed model grows in two, and the second one is the interesting one. A vector of three seniority levels and a vector of two hundred are the same four equations - what changes is how many elements those equations hold, and every sub-element becomes an entity of its own: headcount[l0], headcount[l1], and so on. A model that reads as five lines can be four thousand entities.
This page measures both axes on a workforce aging chain. The scalar benchmark measures the first axis on an ordinary model; read how both are measured first if you intend to compare the two.
The model
A chain of seniority levels. People are hired into it, promoted along it and lost from every level, each level costs its own salary, and three aggregations turn the vector back into single numbers. levels is the parameter the benchmark varies:
from BPTK_Py import Model
def build_chain(levels, stoptime, dt=1.0):
"""A workforce aging chain over `levels` seniority levels."""
names = [f"l{i}" for i in range(levels)]
model = Model(starttime=0.0, stoptime=stoptime, dt=dt, name="chain")
headcount = model.stock("headcount")
headcount.setup_named_vector({name: 100.0 for name in names})
hiring = model.flow("hiring")
hiring.setup_named_vector({name: 0.0 for name in names})
promotion = model.flow("promotion")
promotion.setup_named_vector({name: 0.0 for name in names})
attrition = model.flow("attrition")
attrition.setup_named_vector({name: 0.0 for name in names})
cost = model.converter("cost")
cost.setup_named_vector({name: 0.0 for name in names})
hiring_rate = model.constant("hiring_rate")
hiring_rate.setup_named_vector(
{name: (12.0 if i == 0 else 1.0) for i, name in enumerate(names)})
promotion_rate = model.constant("promotion_rate")
promotion_rate.setup_named_vector(
{name: (0.0 if i == levels - 1 else 0.12) for i, name in enumerate(names)})
attrition_rate = model.constant("attrition_rate")
attrition_rate.setup_named_vector({name: 0.08 for name in names})
salary = model.constant("salary")
salary.setup_named_vector({name: 50000.0 + 1000.0 * i for i, name in enumerate(names)})
# Four equations, written once and holding at every level
hiring.equation = hiring_rate
promotion.equation = headcount * promotion_rate
attrition.equation = headcount * attrition_rate
cost.equation = headcount * salary
# The chain itself is the one part that is written level by level
for i, name in enumerate(names):
arriving = hiring[name] if i == 0 else hiring[name] + promotion[names[i - 1]]
headcount[name].equation = arriving - promotion[name] - attrition[name]
total_headcount = model.converter("total_headcount")
total_headcount.equation = headcount.arr_sum()
total_cost = model.converter("total_cost")
total_cost.equation = cost.arr_sum()
average_salary = model.converter("average_salary")
average_salary.equation = salary.arr_mean()
return model, namesTiming both engines
Both runs go through the ordinary API; only the backend argument differs. Every level’s headcount is requested, so the result frame grows with the model - as it would for a user who actually wants those numbers:
import statistics
import time
import BPTK_Py
def measure(levels, stoptime, runs=3):
"""Median wall-clock time of one full run on each engine, and their deviation."""
equations = [f"headcount[l{i}]" for i in range(levels)] + [
"total_headcount", "total_cost", "average_salary"]
timings, frames = {}, {}
for backend in ("python", "rust"):
samples = []
for _ in range(runs):
model, _ = build_chain(levels, stoptime) # fresh, so no cache carries over
bptk = BPTK_Py.bptk()
bptk.register_scenario_manager({"bench": {"model": model}})
bptk.register_scenarios(scenarios={"base": {}}, scenario_manager="bench")
start = time.perf_counter()
frames[backend] = bptk.plot_scenarios(
scenario_managers=["bench"], scenarios=["base"],
equations=equations, return_df=True, backend=backend)
samples.append((time.perf_counter() - start) * 1000)
timings[backend] = statistics.median(samples)
deviation = max(abs(frames["python"][column] - frames["rust"][column]).max()
for column in frames["python"].columns)
return timings, deviationResults: the time axis
Twelve levels - 279 entities - over a growing number of timesteps. Median of three runs on an Apple M3, Python 3.13:
| Timesteps | Python engine | Rust engine | Speedup | Time saved |
|---|---|---|---|---|
| 400 | 182.2 ms | 4.6 ms | 39.3× | 0.2 s |
| 4,000 | 1,858.6 ms | 31.4 ms | 59.2× | 1.8 s |
| 40,000 | 19,096.5 ms | 303.2 ms | 63.0× | 18.8 s |
Results: the width axis
400 timesteps throughout; only the vector grows. The entity count is what the engine actually loads - five arrayed elements per level, plus the scalars:
| Levels | Entities | Python engine | Rust engine | Speedup | Time saved |
|---|---|---|---|---|---|
| 3 | 72 | 52.7 ms | 2.2 ms | 24.5× | 0.05 s |
| 12 | 279 | 191.5 ms | 5.1 ms | 37.6× | 0.2 s |
| 50 | 1,153 | 779.9 ms | 15.3 ms | 50.9× | 0.8 s |
| 200 | 4,603 | 3,111.6 ms | 61.0 ms | 51.0× | 3.1 s |
The two engines agree exactly. The largest absolute difference across every equation of every row above, both tables, is zero.
What the aggregations cost
arr_sum over two hundred sub-elements is one variadic call in the engine and an operator tree in Python. Running the same models with and without the three aggregations isolates that difference:
| Levels | Python engine | Rust engine | |
|---|---|---|---|
| 50 | with aggregations | 783.2 ms | 15.6 ms |
| 50 | without | 545.5 ms | 14.5 ms |
| 200 | with aggregations | 3,112.4 ms | 61.9 ms |
| 200 | without | 2,185.9 ms | 58.0 ms |
Three aggregations cost the Python engine about 30 % of the whole run, at both widths. They cost the Rust engine 6-7 %.
Reading the tables
An arrayed model gains far more than a scalar one. The scalar benchmark measures 6-10× on the same machine; this one measures 24-63×. The reason is what an array is on the Python side: every sub-element carries its own operator tree and its own memoised lambda, so a vector of 200 is 200 times the interpreter overhead. The engine sees 200 flat entities and evaluates them in a loop.
The ratio grows with size instead of shrinking. In the scalar benchmark the speedup falls at the very top, because converting the result to a DataFrame costs the same on both paths and eventually dominates. Here the evaluation is so much heavier on the Python side that it keeps dominating: the ratio climbs from 24× to about 60× and then settles. What settles it is the same DataFrame conversion - at 200 levels the frame has 203 columns.
Width costs what timesteps cost. 3 levels to 200 is a factor of 64 in entities and a factor of 59 in Python’s runtime, near enough proportional; the Rust engine grows by a factor of 28 over the same range, because a fixed setup cost is a larger share of its two milliseconds. There is no width at which arrays become disproportionately expensive - they are ordinary entities to both engines.
Serialising the model stays negligible. Turning a 4,603-entity model into JSON takes 3.3 ms, and 279 entities take 0.3 ms - against a run of seconds. This was the one place where arrays might plausibly have cost something the scalar benchmark does not see, and they do not.
Where this leaves the advice. For a scalar model the Rust engine pays off with size and repetition. For an arrayed one it pays off almost immediately: at twelve levels and 400 timesteps - a small model by any measure - the run is already 39× faster, and the absolute saving becomes seconds as soon as either axis grows.