The same declared model can be run by five solvers ('ode', 'binomial', 'gillespie', 'tau', and 'abm'), some with options. This page explains how Starsim chooses between them, the engines and kernels that do the work, the tests that pin each path, and how the engines compute transmission and apply the routes, effects, interventions, and results that are kept separate from the disease. For the solvers from a modeler’s side (what each is good for, and how their results differ), see Solvers and backends.
How the solver is chosen
The default, solver='auto', picks the simplest, fastest representation that runs the model without changing what it means: an ODE if every part of the model has a compartmental form, and agents otherwise. The line the sim prints when it is initialized says which it chose, and sim.solver_reason (and describe()) says why:
import starsim as ssdefmake_sir(dur=ss.days(10)): S, I, R = ss.states('S I R')return ss.Disease(infection= S >> I @ ss.perday(0.3) * I,recovery= ss.Transition(I, R, dur=dur),init_prev=0.01, )kw =dict(unit='day', stop=100)sim = ss.Sim(make_sir(), **kw).init() # Homogeneous mixingprint()sim = ss.Sim(make_sir(), networks=ss.RandomNet(n_contacts=10), **kw).init() # A network of individualsprint()sim = ss.Sim(make_sir(ss.lognormal(mean=ss.days(10), std=ss.days(4))), **kw).init() # A lognormal dwell time
Sim: sir, 10,000 people, as an ODE, days 0 to 100
sir: R0 = 3, with 100 infectious at the start
Sim: sir, as 10,000 agents, days 0 to 100
Running as agents: RandomNet is an individual-level network
Note: consider a shorter step (e.g. dt=ss.days(0.25)): each step is 1 day, and infection transmits at up to 0.3 per day (per infectious contact), 26% in a step. Taking whole steps, agents may differ from the ODE by an amount of that order; sim.check_dt() shows how much it matters
Sim: sir, 10,000 people, as an ODE, days 0 to 100
sir: R0 = 3, with 100 infectious at the start
Note: recovery has a lognormal dwell time (nearly fixed, not exponential), approximated by 6 Erlang stages (mean 10 days kept; std 4.08 days, against 4 days given)
The rule behind this is that 'auto' never approximates the mixing and never drops anything individual-level: a model runs as agents if it has an individual-level network, n_agents=, a function of uids, a module written as code that keeps a value for each person, or anything that needs individual histories. Being written as code isn’t a reason by itself: a module that keeps nothing per person runs beside the compartments. The full list is in What needs agents, and the rules themselves are in solvers/choose.py (below).
The one thing 'auto' approximates unasked is a dwell time that isn’t exponential, Erlang or fixed, as in the third example. Sending a simple model to agents over one lognormal would be the bigger surprise, so it uses Erlang stages that keep the mean exactly, and prints the standard deviation they give. A dwell time too variable for stages (a coefficient of variation above 3) makes 'auto' choose agents.
sim.describe() lists each ingredient that matters for the choice and whether it runs as compartments: yes, approx (and how), or no (and why). Code is never analyzed in advance: a function is called with what the solver has (one value per agent, or one per group), and if it asks for something only agents have, the run stops with an ss.CapabilityError that names it.
'auto' only ever chooses ode or abm; a stochastic compartmental solver, exact dwell times (ss.solvers.ode(dwells='exact')), or event-driven agents (ss.solvers.abm(step='events')) are used only when asked for. Any solver can also run the model as a branching process (branching=True; see Branching and rescaling). Every solver and its options are listed in Solvers and backends.
Approximations are announced
When a solver is asked for that the model doesn’t natively fit, Starsim uses the standard approximation where one exists, and says so in one line. It stops only when there is no sensible approximation. For example, a random network run as an ODE becomes homogeneous mixing at the network’s mean number of contacts:
>>> ss.Sim(sir, networks=ss.RandomNet(n_contacts=10), solver='ode', unit='day', stop=100).init()
Sim: sir, 10,000 people, as an ODE, days 0 to 100
sir: R0 = 30, with 100 infectious at the start
Note: RandomNet approximated by homogeneous mixing with 10 contacts per day (its mean number of contacts)
Each ingredient has one standard approximation (lasting partnerships use the Anderson–May partner-change model; a distribution on beta uses its mean). Things with no compartmental form, such as contact tracing, stop with an ss.CapabilityError before the run, saying what to do instead.
By default each approximation prints a note; approximations='allow' keeps the notes without printing them, and approximations='error' makes any approximation an error (under 'auto', that sends a lognormal dwell time to agents). Either way the notes are kept in sim.notes, so they survive into saved output.
What decides which engine runs a model
One ordered list of rules decides it, in solvers/choose.py, read for solver='auto', for the CapabilityError of a solver that was asked for, and for why a run is stepped (in sim.describe()). Each rule says what it looks for, the engines that can run a model with it, and the message; a module adds its own reasons with agent_reasons(), as the built-in ones do. In the table, compartments is one step at a time, kernels is the whole run in one call, and jax is the ODE in JAX.
ss.solvers.solver_table()
what
agents
compartments
kernels
jax
0
n_agents was given
yes
if asked
no
no
1
people were given
yes
if asked
no
no
2
How often or how long ago an event happened, w...
yes
no
no
no
3
A step triggered by a product delivered by sev...
yes
no
no
no
4
What a module says it needs agents for: see `a...
yes
no
no
no
...
...
...
...
...
...
58
solver="..." is stochastic (gradients are for ...
yes
yes
yes
no
59
Whatever stops the kernels taking the whole run
yes
yes
yes
no
60
"..." has a movement register
yes
yes
yes
no
61
A probability of 1 over a period at the start,...
yes
yes
yes
no
62
the rate of the transfer "..." of "..." change...
yes
yes
yes
no
63 rows × 5 columns
The engines
Each declared module gets one engine when the sim is initialized. A module written as code, with no declared states or transitions, gets an empty engine and runs its own step(). Beside the engines, solvers/deliveries.py is what an intervention does differently as compartments (it reaches a fraction of each group, not a list of agents), and solvers/particles.py steps many runs of one program together for the particle filter.
Expands states into compartments (one per group, branch and Erlang stage) and transitions into a flat list of flows, the program. The kernels (small functions in solvers/backends/kernels.py, compiled with Numba, or their Rust copies) step that program for all the compartmental solvers. mixing.py turns networks and mixing pools into contact rates between groups.
Each state is a boolean array over agents. A dwell time is drawn on entry and the exit scheduled; competing rates out of a state are resolved in one draw; transmission is drawn over each network’s edges. The moves of a step are taken in rounds (below), each a few NumPy passes over the agents who move. The heavy passes over agents and edges are compiled.
JointStates
solvers/joint_states.py
Diseases tracked jointly, when a condition in one refers to a state of another (ss.modify(hiv.rel_trans, multiply=3, who=syph.I)). Their compartments are crossed and run as one program, up to a budget (ss.options.max_joint); above it, 'auto' uses agents.
Three details of the compartmental engine explain most of its behavior:
A whole run in one call. If nothing needs Python between steps (no module written as code with a step of its own, such as an intervention, analyzer, or connector, and no modifier with who=, when=, or a function as its amount), the kernels run the whole simulation in one call; otherwise they are called once per step. Values that vary in time (ss.timeseries, ss.seasonal, and ss.by of them, by place or age) are tabulated for every step before the run, and so are deliveries on a schedule to people chosen by state and group (ss.Delivery.tabulated), including a move between states such as a seed (ss.Event(S >> I, when=...)), so none of these makes a run stepped. The results are the same, and sim.describe() says why a run is stepped.
Rust, Numba, or plain. The kernels run in Rust, compiled by Numba, or as plain Python, with the same results (see the kernels, below).
JAX for gradients and fitting.internals/gradients.py (ss.gradients) runs the same program with the same steps in JAX, so results can be differentiated (sim.gradient()) and many parameter sets run at once. ss.solvers.choose_engine() uses JAX for a run, or for a sweep of many parameter sets in one call (ss.parallel(sim, pars=...)), only with sim.jax = True; gradients and sim.fit() use it on their own. sim.engine_reason says which ran, once the sim has been run.
What a step does, under each solver (competing exits, when rates are read, cascades, and the order between modules), is in From declaration to run; this is how the agent engine does it. It takes the transitions that are not transmission first, then transmission, each as one sweep over the disease’s transitions, in an order worked out once from the graph (Agents.sweep()): each transition after those that bring people into its source, a trigger (when=) after the event it waits for, and the triggers out of a state ahead of its scheduled exits; a loop of states is broken at the transition declared first. At each transition in turn, those it moves now leave (by its trigger, or because their exit by it is due, including people who arrived in its source earlier in the sweep with a wait that rounds to no steps), and everyone who left by it is sent to their destinations (a transition with several destinations takes one draw for each, and ss.stay puts people back); each arrival is dated, and has its branch, dwell time and exit drawn if its state schedules them. Someone who reaches a transition the sweep has already passed is taken by another sweep, until nobody moves, so the moves with no wait happen in the step that caused them, as in the ODE, and in the order a modeler writing the code by hand would take them (Agents.finish_moves()). The step’s infections are dated at its end, and start a cascade of their own. Every way in from code (an intervention’s move, ss.move(), disease.infect()) is a round of its own. The compiled passes take the same rounds.
Each of the engine’s draws is keyed by its stream, the step, and the round (Agents.draw(), Agents.new_round()): each transition of the sweep is a round, as is each way in from code. A draw is not keyed by how many times the stream has been called, so it does not depend on how the work of a round is batched; an agent moves at most once in a round, so each agent’s draws from a stream are its own (STARSIM_CHECK_DRAWS=1 checks it).
The random draws are shared too: solvers/backends/draws.py has exact draws from counter-based uniforms, so a stream doesn’t depend on how work is split across threads. Users’ own code can be compiled with @ss.jit (solvers/backends/jit.py), with the same results as when run as Python.
The kernels
The same kernels exist as several sets, which give the same numbers: Rust’s (solvers/rust/, swapped in by solvers/backends/rust.py), Numba’s (solvers/backends/kernels.py, passes.py, events.py, draws.py), and plain Python’s (the compartmental kernels run without compiling; for the ODE, solvers/backends/ode.py in NumPy). ss.options.kernels chooses between them; by default Rust is used where it is built, then Numba, with plain Python for compartmental runs too short to be worth compiling. JAX is not one of these: it is an engine for a whole run of the ODE (sim.jax). Numba’s kernels are the reference: a change is made there first, then in Rust’s (solvers/rust/README.md says how), and test_kernel_sets.py checks that every set still gives the same numbers. What each set works with is in Solvers and backends.
Execution paths and the tests that pin them
Each path is checked against a reference where one exists: the ODE against SciPy and known answers, and the other paths against the ODE or against a path already checked. The developer guide lists each path, when it is used, and the tests that pin it.
How the engines compute transmission
A bare state in a rate means contact with people in that state, through whatever routes the sim has (all passed with networks=). With no network, beta * I is the textbook beta * I / N, and beta is a transmission rate; with a network or mixing pool, beta is a rate per contact. It is one rule: the hazard is beta times the contact rate times the fraction of contacts who are infectious, and the default contact rate is 1. I.n is used literally, except that a contact matrix acts on it in an acquisition’s rate; this is how environmental and vector-borne transmission are written, with the reservoir a state declared with people=False. A list of routes is a layered population, whose forces of infection sum, and space is a stratification by patch, coupled by a mobility model. For the modeler’s view, see Mixing.
Each engine computes the rule in its own terms, from the same declaration:
As compartments, the routes are turned into contact rates between groups when the sim is initialized (ss.solvers.Mixing, in solvers/mixing.py): homogeneous mixing is one contact per unit of time, a mixing pool sets the number of contacts, mixing pools between age groups give a contact matrix, and space adds the movement between patches. The force of infection on each group is then beta times the contact rates times the fraction infectious in each group it meets, which the kernels evaluate as part of the program. Nothing is drawn per contact.
As agents, transmission is drawn over each network’s edges (ss.solvers.Transmission, in solvers/agents.py): for each edge with an infectious person at one end and a susceptible at the other, the chance in the step is 1 − exp(−beta·dt·rel_trans·rel_sus), with the effective rel_trans of the source and rel_sus of the target. Every factor of a contact’s hazard is in the exponent: beta, the source state’s coefficient, rel_trans, rel_sus, the edge’s weight, and the modifiers on the route. So a contact infects at most once in a step, and an effect gives the same answer wherever it is written (halving rel_sus is the same as halving beta). One consequence: a chance of 1 per contact stays 1 under any factor above 0, since it is an infinite hazard. A network with acts per partnership works out a hazard per act, and the acts in the step scale it the same way. The routes are taken in the order of sim.networks, and someone reached by several sources is infected by the first, who is recorded as the source in the infection log. A mixing pool as agents is a force of infection on each agent from the fraction infectious in the pool, with no edges.
A disease written as code (ss.AgentDisease) transmits through the same ss.solvers.Transmission code, so it sees the same modifiers and products, gives the same edges the same random numbers, and attributes infections the same way; a bare beta is a rate in the sim’s unit, as in a declared disease, and on a network with acts a probability per act.
Each edge’s random number is drawn by its source and target, from a stream named after the disease, so adding a network or another disease doesn’t change which edges transmit (see Random numbers).
Effects, interventions, and results
Each of these is a separate object from the disease, so changing one doesn’t touch it. A few rules hold across them; for how to use each, see Connectors, Interventions, Products, and Results.
Effects. One module changes another either with ss.modify (while a condition holds, set, multiply, or add to a named quantity) or with ordinary Python in its step(). A modifier never overwrites what it changes: each modifiable quantity has a base value (what its module set) and an effective value (the base with every active modifier applied). Reading or assigning sir.rel_sus[uids] uses the base, and sir.current_value('rel_sus', uids) is what transmission will use, through ss.solvers.Transmission, shared by declared diseases and those written as code. When a modifier switches off, the quantity is back at its base, with nothing to restore, and several modifiers on one quantity combine by one fixed rule, whatever their order: set, then multiply, then add. Products are modifiers too, so sim.what_modifies(sir.rel_sus) lists vaccines with the rest, before anything runs. Conditions (who=, e.g. ss.person.age > 65, an ss.conditions.Cond) are lazy, and tried on whatever the solver provides: age > 65 works on groups if 65 is a bin edge, and otherwise needs agents. A modifier that no solver would apply raises a ValueError when the sim is built. Code that writes another module’s arrays directly is allowed, but isn’t tracked.
Interventions. A product (ss.Vaccine) is what is given, and a delivery (ss.vaccinate) is who gets it, when, and how many. A delivery’s outcome such as screening.positive is an event, and a cascade of care is a list of deliveries each triggered by another’s event. Each delivery runs after the deliveries whose events trigger it, whatever order they are listed in, so a follow-up acts in the same step as its trigger. As compartments a delivery reaches a fraction of each group (solvers/deliveries.py), with tests as flows thinned by sensitivity and specificity; a delay, the joint of two events, .since, or trace= need agents.
Results. Every state gets a count and every transition a flow of the same name. Each result knows its kind (flow, stock, or ratio), which sets how it is summarized, resampled, and compared with data, so no caller needs to know which it is; a ratio is never averaged across groups. An ss.MultiSim’s results are the same ss.Result objects with a row for each run, so the same methods work on one run or many.
The code is in starsim/networks.py and solvers/mixing.py (routes), connectors.py and internals/conditions.py (effects), interventions.py, products.py and solvers/deliveries.py (interventions), and results.py.
Limits
Solvers that take whole steps differ from the ODE by an amount of order dt, which is why the agreement tests use a quarter-day step. With whole weekly steps (ss.solvers.binomial(substeps=1)) on an SIR with daily rates, the chain binomial’s peak is about a third higher than the ODE’s; by default it takes enough substeps to stay within a few percent. A compartment holds about x/2 too many people, for x its total rate out times the substep, so Compartmental.binomial_substeps() keeps x to 5%, and near the epidemic threshold (where a small bias in R moves prevalence a lot) to a tenth of |R0 − 1|.
Agreement is tested on simple models. For complex ones (a sexual network as an ODE), solvers are expected to differ; the guarantee is that Starsim says what it approximated.
Gillespie is exact but takes events one at a time; for millions of events in groups that mix, 'tau' is the practical choice.
Gradients exist only for solver='ode' and for rates written as expressions. One model can’t yet mix solvers (an ODE with rare importations as discrete events).
Modifiers on one quantity combine by the one rule in Effects (set, then multiply, then add). Protection from several strains a person holds needs a different rule (the strongest, or a capped sum), which is ss.immunity; anything else is step() code.
Functions in who= and ss.modify keep no memory apart from until=; anything that must remember is a class.
Routes and Effects: the design record’s account of transmission and modifiers
dust2: Introduction: the engine behind odin2, which runs one compiled model deterministically or as many stochastic particles, as Starsim’s compartmental engine does
GillesPy2: exact (Gillespie) and tau-leaping solvers for stochastic compartmental models, with several compiled backends