Benchmark: A Scalar Model

How much the Rust engine is worth on a scalar System Dynamics model, measured from 400 to 400,000 timesteps.

Keywords

benchmark, performance, rust, execution engine, system dynamics, bptk, bptk-py, python

Benchmark: A Scalar Model

What the choice of engine is worth on an ordinary, scalar model, with a model small enough to print in full — copy the four blocks into a file and you have the benchmark. The arrayed benchmark measures the other axis.

The model

An SIR epidemic with vaccination and a hospital-capacity feedback. People move from susceptible to infected, from infected to recovered or deceased, and from susceptible to vaccinated; the recovered slowly lose their immunity and return. Mortality is not a constant: it rises once the infected outnumber the hospital beds. That feedback, plus the waning immunity, is what makes the model oscillate rather than settle after one wave — and what makes it worth simulating rather than solving.

Five stocks hold the population:

from BPTK_Py import Model
from BPTK_Py.sddsl import functions as sd


def build_sir_model(stoptime=100, dt=0.25):
    model = Model(starttime=0, stoptime=stoptime, dt=dt, name="sir")

    susceptible = model.stock("susceptible")
    infected = model.stock("infected")
    recovered = model.stock("recovered")
    vaccinated = model.stock("vaccinated")
    deceased = model.stock("deceased")

    susceptible.initial_value = 9990.0
    infected.initial_value = 10.0
    recovered.initial_value = 0.0
    vaccinated.initial_value = 0.0
    deceased.initial_value = 0.0

Five flows move people between them. Each stock equation is simply the flows in minus the flows out — the arithmetic operators build the equation graph, they do not compute anything yet:

    infection = model.flow("infection")
    recovery = model.flow("recovery")
    vaccination = model.flow("vaccination")
    death = model.flow("death")
    waning_immunity = model.flow("waning_immunity")

    susceptible.equation = -infection - vaccination + waning_immunity
    infected.equation = infection - recovery - death
    recovered.equation = recovery - waning_immunity
    vaccinated.equation = vaccination
    deceased.equation = death

Seven constants parameterise it:

    contact_rate = model.constant("contact_rate")
    transmission_prob = model.constant("transmission_prob")
    recovery_time = model.constant("recovery_time")
    mortality_rate = model.constant("mortality_rate")
    vaccination_rate = model.constant("vaccination_rate")
    immunity_duration = model.constant("immunity_duration")
    hospital_capacity = model.constant("hospital_capacity")

    contact_rate.equation = 8.0
    transmission_prob.equation = 0.03
    recovery_time.equation = 14.0
    mortality_rate.equation = 0.001
    vaccination_rate.equation = 0.01
    immunity_duration.equation = 180.0
    hospital_capacity.equation = 200.0

And three converters carry the parts that are neither a level nor a rate. capacity_pressure is the feedback: at or below capacity it is 1 and mortality is the base rate; above it, it grows with the number of infected and takes effective_mortality with it.

    total_population = model.converter("total_population")
    capacity_pressure = model.converter("capacity_pressure")
    effective_mortality = model.converter("effective_mortality")

    total_population.equation = susceptible + infected + recovered + vaccinated

    infection.equation = sd.max(
        0,
        contact_rate * transmission_prob * susceptible * infected / total_population,
    )
    recovery.equation = sd.max(0, infected / recovery_time)
    vaccination.equation = sd.If(
        sd.time() > 10,
        sd.min(susceptible * vaccination_rate, susceptible / model.dt),
        0,
    )

    capacity_pressure.equation = sd.max(1, infected / hospital_capacity)
    effective_mortality.equation = mortality_rate * capacity_pressure
    death.equation = sd.max(0, infected * effective_mortality)
    waning_immunity.equation = sd.max(0, recovered / immunity_duration)

    return model

Twenty elements in total, thirteen of them tracked in the runs below. It is deliberately a small model: what changes from row to row is only how many timesteps it is run for, from 400 to 400,000, because that is what a parameter sweep, a Monte Carlo run or a long horizon actually does to you.

Timing both engines

Both runs go through the ordinary API. Only the backend argument differs:

import time
import BPTK_Py

EQUATIONS = [
    "susceptible", "infected", "recovered", "vaccinated", "deceased",
    "infection", "recovery", "vaccination", "death", "waning_immunity",
    "total_population", "capacity_pressure", "effective_mortality",
]


def measure(stoptime, dt=0.25, runs=3):
    """Median wall-clock time of one full run, on each engine."""
    model = build_sir_model(stoptime=stoptime, dt=dt)

    bptk = BPTK_Py.bptk()
    bptk.register_scenario_manager({"bench": {"model": model}})
    bptk.register_scenarios(scenarios={"base": {}}, scenario_manager="bench")

    times = {"python": [], "rust": []}
    results = {}
    for _ in range(runs):
        for backend in ("python", "rust"):
            model.reset_cache()   # without this the second run reads the first one's memo
            start = time.perf_counter()
            results[backend] = bptk.run_scenarios(
                scenario_managers=["bench"], scenarios=["base"],
                equations=EQUATIONS, backend=backend,
            )
            times[backend].append(time.perf_counter() - start)

    median = {b: sorted(v)[runs // 2] for b, v in times.items()}
    largest_difference = (results["python"] - results["rust"]).abs().max().max()
    return median["python"], median["rust"], largest_difference


for stoptime, steps in [(100, "400"), (1_000, "4,000"),
                        (10_000, "40,000"), (100_000, "400,000")]:
    python_s, rust_s, difference = measure(stoptime)
    print(f"{steps:>9} steps   python {python_s * 1000:9.1f} ms"
          f"   rust {rust_s * 1000:8.1f} ms   {python_s / rust_s:5.1f}x"
          f"   largest difference {difference:.1e}")

reset_cache() is not decoration. Without it the second engine reads the values the first one memoised and reports a speed that belongs to neither.

Results

Median of three runs on an Apple M3, Python 3.13:

Timesteps	Python engine	Rust engine	Speedup	Time saved
400	26.0 ms	3.9 ms	6.7×	22 ms
4,000	257.5 ms	24.3 ms	10.6×	0.2 s
40,000	2,502.5 ms	260.5 ms	9.6×	2.2 s
400,000	27,548.4 ms	4,398.0 ms	6.3×	23.2 s

The two engines agree exactly: the largest absolute difference across all thirteen equations, in every row, is zero.

Your own numbers will differ — the ratio is what to compare, not the milliseconds. One caveat if you build the engine from source: maturin develop produces a debug build, roughly five times slower than the release build that ships in the wheel. Benchmark against maturin develop --release, or against an installed pip install BPTK-Py.

Reading the table

The ratio is not the interesting column; the last one is. At 400 timesteps the Rust engine saves 22 milliseconds, which nobody notices. At 400,000 it saves 23 seconds — the difference between a run you wait for and a run you switch away from. Put that inside a hundred-iteration sweep and it is the difference between forty minutes and six.

The ratio drops at the very top, and the engine is not the reason. Simulating 400,000 steps costs the Rust engine about 1.9 s; turning the result into a pandas DataFrame costs about 2.6 s — and that conversion is the same Python code on both paths. The simulation itself stays 15–20× faster throughout; what you measure end-to-end is diluted by the part neither engine can speed up. Serialising the model to JSON and parsing it in Rust, by contrast, cost about 0.1 ms each and never matter at any size.

Which is why the advice under Troubleshooting is what it is: the Rust engine pays off with size and with repetition. A small model run once is dominated by setting it up.