From declaration to run

🧠0% 🤖100%

A disease written as code, as a class whose step() method changes arrays of agents, can be run but not read: Starsim can’t describe the model, draw it, or run it any other way. A declared disease is instead a set of objects (states, transitions, parameters) that describe the model without running it. Everything else is derived from them: the description, the diagram and equations, the stratified model, the exports to other frameworks, and the engines that run it as an ODE, a stochastic compartmental model, or agents.

This page follows a model from those objects to a run: the objects, how they are resolved and checked, what sim.init() does, how time and the order of a step work, and how Starsim reports what it assumed. For how to write a disease as a modeler, see Diseases and Writing models.

Declaring a model

States and transitions

A state (ss.State) is something a person can be in: a compartment in an ODE, a boolean array over agents in an agent-based run. A transition (ss.Transition) moves people between states, driven by a rate, a dwell time, a branch probability, or a trigger. The ways of writing a model (operators, ss.Transition(), a class body, or define_states() and define_transitions() in __init__(); see Writing models) all build the same objects, and give the same results. This is a rule of the design: a keyword means the same thing in every form, and subclassing adds behavior, not vocabulary. A transition takes the state objects themselves, never their names as strings, so nothing is inferred from a name. What is declared in __init__() is collected as it is given and built once, when the outermost __init__() returns (BuiltAfterInit in starsim/diseases.py), so a subclass can add to its parent’s states before anything is built; after that, declaring more raises an error.

Rates are expressions

When you write 0.3 * I, nothing is computed. I is a state, so the result is an expression (ss.expressions.Expr, in starsim/internals/expressions.py): a small tree whose leaves are parameters, states, counts (I.n), and properties of people (ss.person.age), joined by arithmetic and NumPy functions. Each engine evaluates the same tree in its own terms (a value for each group, for each agent, or as JAX arrays for gradients), which is what lets one model run under any solver.

Two kinds of leaf differ. I.n is the number of people in I, used literally (or, in an acquisition’s rate with a contact matrix, the number met through it). A bare I means contact with people in I: with no network that’s I.n / N, so beta * I is the textbook beta * S * I / N, and with a network each engine works it out from who actually meets whom, which no single formula can say. Since that is a sum over contacts, a rate must be linear in its bare states (beta * I * I raises a ValueError that suggests the two ways to say what was probably meant). Engines and transmission covers how each engine turns a bare state into a force of infection.

The kind of a transition is derived from its rate, not declared: if the rate mentions a bare state, the transition is contact-mediated and that state is a source of infection; otherwise it is spontaneous. For sir = ss.Disease(infection=S >> I @ 0.3 * I, recovery=I >> R @ ss.days(10)):

sir.infection.contact     # True: its rate mentions I
sir.infection.terms       # {I: par(0.3)}: the coefficient of each source of infection
sir.recovery.contact      # False

Parameters

A number written inline becomes a parameter of the module, named after its slot (beta for the coefficient of I, dur_I for the time in I; the rules are in the appendix of Writing models). The slot decides what a value means, and the form is free: a value may be a number, an ss.par() with bounds or a source, or vary by group, in time, or between people. See Parameters.

Dimensions

Age, sex, place, and risk group are all handled by one operation, stratify(), which divides a finished model by a property of people. It doesn’t copy the states: the dimension becomes an axis, with one row per group in a compartmental run, and an array on the people as agents. The base model and the stratification don’t refer to each other, so one disease can be reused with different structures, and a stratified result is still one result (res.I.by('age')). The code is in starsim/internals/dimensions.py (ss.dimensions.Strat), and the groups a sim ends up with are sim.groups. See Stratification.

The vocabulary

The design record lists every public name in one place (The complete vocabulary), and a name is added only after asking whether Starsim already has one, whether it can be written with what exists (a bednet is an ss.modify), whether more than one kind of model would use it, and whether a keyword is enough. Modules, containers, and base classes are CamelCase (ss.Sim, ss.Vaccine); the pieces put inside them are lowercase (ss.days, ss.modify, ss.vaccinate). Importers write only names a person would write by hand.

The graph: resolve, complete, check

When ss.Disease(...) is called, it builds the module’s graph (ss.graph.Graph, in starsim/internals/graph.py). The graph is where the objects you wrote become a model:

Collect and name. It gathers the states, transitions, parameters, and results, and names any that have no name: from the keyword they were passed under, then a name= argument, then the variable they’re bound to in the calling code. A transition with none of these is named infection if it is the disease’s only transmission (as in ss.SIR), and otherwise gets a name built from its endpoints (I_to_R); a state with none raises a ValueError, since a made-up state name would leak into results. Names are picked up from variables in this one place; it can miss (a transition written inline in the call gets the fallback name), and describe() shows the names used.
Resolve. It finds which states are people and which are plain counts (a water reservoir, a mosquito population), derives the kind of each transition from its rate, marks the sources of infection, works out which transitions are acquisitions of the disease, and finds the first state (where people start and are born into).
Complete. It turns inline values into named parameters, and groups the exits of each state into those that compete (rates, resolved together in one draw; see below) and those that split on entry (p=).
Check. It checks what can be checked without a sim: the time dimension of every rate, that branch probabilities add up, that exposure is linear, and that keywords name real parameters.

So mistakes in the model’s structure are errors when the disease is built, at the line that made it. What needs the sim (the unit of bare numbers, values against their bounds, the shapes of values by group, the choice of solver) is checked at sim.init().

Guessed defaults

Short models leave a lot unsaid: how a dwell time is distributed, how people mix, how many there are. Starsim fills each gap with a guessed default, and describe() marks every guess:

import starsim as ss

S, I, R = ss.states('S I R')
beta = ss.par(bounds=[0.1, 1.0], source='to be fitted')
sir = ss.Disease(
    infection = S >> I @ beta * I,
    recovery  = I >> R @ ss.lognormal(mean=ss.days(10), std=ss.days(3)),
    init_prev = 0.01,
)
sim = ss.Sim(sir, stop=100)
sim.describe()

Sim 1 module, 10,000 people [guessed], days 0 to 100

solver	ode	guessed: every part of the model has a compartmental form
unit	day	from the typed values in the model

Disease 'sir' 3 states, 2 transitions

states	S, I, R
	I infectious
	init: init_prev = 0.01; the rest in S

infection	S -> I transmission, from I
	beta = 0.316228 per day, between 0.1 and 1	guessed: the middle of its bounds
	source: to be fitted
	rate = 0.316228 × I.n / N	homogeneous mixing
recovery	I -> R spontaneous
	time in I ~ lognormal(mean=10 days, std=3 days), as 11 Erlang stages	guessed: stages from the mean and standard deviation

Mixing and interventions

mixing	homogeneous, 1 contact per day	guessed: no network given

Results and running

results	S, I, R, and one flow per transition; also infected, prevalence
plan	sir: as compiled code, the whole run in one call
schedule	every day: sir

Notes

note	unit='day' (from the typed values in the model)
note	recovery has a lognormal dwell time (nearly fixed, not exponential), approximated by 11 Erlang stages (mean 10 days kept; std 3.02 days, against 3 days given)

Per step of sir 1 day

infection	beta = 0.316228 per day → 0.316228 per step
recovery	rate 1.1 per day in each of 11 stages → p = 1 - exp(-1.1 × 1) = 0.6671 per step

Every [guessed: ...] is a choice you didn’t make: 10,000 people, the ODE solver, homogeneous mixing, 11 Erlang stages for the lognormal dwell time, and beta starting at the geometric middle of its bounds. Each guess is tested both ways: a model that makes it has it marked, and one that doesn’t, doesn’t.

Because the model is data, its other views (equations(), diagram(), and the table of transitions, all in starsim/internals/display.py; see Writing models) and its exports (below) are read from the same graph.

What `sim.init()` does

Creating a sim only stores what it was given (as copies, so one sir can go into several sims). The rest happens at sim.init(), which sim.run() and sim.describe() call if needed. In order (Sim.init() in starsim/sim.py, and choose_solver() and the steps that build a declared model in starsim/solvers/choose.py):

Values. Values given by name (pars=dict(beta=0.1)) are set, and each parameter is checked against its bounds (a ValueError).
Groups. The stratifications of the sim and its modules are collected into sim.groups.
Unit. The unit of bare numbers is chosen (below).
Solver. Starsim lays the model out as compartments (which runs none of its code) and, under solver='auto', chooses an ODE if every part has a compartmental form, agents otherwise; a compartmental run reuses that layout. See Engines and transmission.
Objects. The timeline, the people, and for each declared disease an engine, the object that runs its transitions under the chosen solver. This is the one place where the solver changes which objects exist.
Starting values. Random number streams (Random numbers), people’s attributes, then each state’s initial conditions.
Results and the loop. Results are allocated, data is attached to them by name, and the loop’s functions are collected. Parameters are locked, so setting one that doesn’t exist raises an error.
Report. One line saying what will run and how, a note for each inference or approximation (kept in sim.notes), and for a compartmental model with transmission, R0 and the number infectious at the start. The reason for the solver is in sim.solver_reason, and in describe().

S, I, R = ss.states('S I R')
sir = ss.Disease(
    infection = S >> I @ 0.3 * I,
    recovery  = ss.Transition(I, R, dur=10),
    init_prev = 0.01,
)
sim = ss.Sim(sir, start='2020-01-01', stop='2020-07-01')
sim.init()

Sim: sir, 10,000 people, as an ODE, 2020-01-01 to 2020-07-01
  sir: R0 = 3, with 100 infectious at the start

Sim(n=10000; 2020-01-01—2020-07-01; diseases=sir)

Errors are raised at the earliest point that can detect them: a mistake in the model’s structure when the module is built, a value out of bounds at step 1, and a function that asks for something the solver can’t give when it is first called. After sim.init() the sim’s settings are frozen: a changed copy is made with sim.new().

The runtime objects

The engines and modules are covered in Architecture and Engines and transmission. The other objects a module written as code reads and writes are:

People. As agents, sim.people (starsim/people.py) has one entry per agent: age, sex, whether alive, and every module’s states. In a compartmental run it has one row per group, so people.age in a model not stratified by age raises an ss.CapabilityError saying to stratify or run as agents.
Arrays. As agents, a declared state is an ss.BoolState over agents (starsim/internals/arrays.py), beside the ss.FloatArr and others that code adds. Arrays grow as agents are born; the dead leave at the end of the step. Writes by code to a declared state (sir.I[uids] = True) are tracked, so counts see them; writes through .raw, the NumPy array underneath, are not. See Array indexing.

Units, time points, and steps

Durations, hazard rates, probabilities per period, and frequencies are types of their own (starsim/internals/time.py, with dates in internals/dates.py): converting one to a timestep is the type’s job (a hazard r becomes the probability 1 − exp(−r·dt), a frequency the expected count f·dt), and arithmetic that mixes them wrongly raises a TypeError. A bare number is typed by its slot (after @, a rate; for dur=, a duration), in the sim’s unit, which is taken from unit= or inferred from the typed values and dates (Time lists the order). The unit is never guessed from the size of a number (0.05 is plausible per day and per year); Time says when inferring it prints a note.

A sim has time points 0, 1, ..., N, from start to stop. Step k runs from time point k to time point k + 1:

time point   0         1         2               N
             |--step 0-|--step 1-|--  ...  --|--step N-|
stocks       I[0]      I[1]      I[2]        I[N]
flows        inf[0]    inf[1]    inf[2]      inf[N]

A step is taken from every time point, the last one included, so every time point is a period:

A stock at time point k (e.g. I[k], the number infectious) is the state at the start of step k. I[0] is the initial condition.
A flow at time point k (e.g. infection[k], the new infections) is the total during step k.

So with yearly steps, the row for 2020 holds the prevalence on 1 January 2020 and the infections during 2020, which is how data labeled “2020” are usually meant. The state after the last step isn’t recorded, so after a run sim.people is one step past the last time point.

Throughout step k, in every module, sim.ti is k (the index, where results go) and sim.now is time point k as a date or a duration (the value, at which values that vary in time are read).

Because a stock is the state at the start of a step and a flow is the total over it, each stock is the one before it plus the flows of the step between them:

import numpy as np

sim = ss.Sim(diseases=ss.SIR(), total_pop=10_000, unit='day', stop=60, solver='binomial', verbose=0)
sim.run()
r = sim.results.sir
print('I[0], the initial condition:', r.I[0])
print('S[k+1] - S[k] == -infection[k] at every step:', np.array_equal(np.diff(r.S), -r.infection[:-1]))

I[0], the initial condition: 100.0
S[k+1] - S[k] == -infection[k] at every step: True

Why stocks and flows are recorded separately

A stock is a snapshot and a flow is a sum over an interval, so they are ready at different times. The stocks are ready at the start of a step, before anything has changed. The flows are ready only at its end, after every module has run and the deaths have happened. One call to record both can get only one of them right. Hence two methods:

record_states() runs at the start of each step, before any module steps.
record_events() runs at the end, after the deaths and before the analyzers.

Starsim 3 had one method, update_results(), called after the step, so its stocks were the state at the end of the step: Starsim 3’s n_infected[k] is Starsim 4’s I[k + 1] (see Migrating from v3). How a module written as code fills each is in Custom modules.

The order of a step

Each step, the loop calls the modules in Starsim’s long-standing order (starsim/internals/loop.py):

Order	Who	What happens
1	sim, then every module	`advance()` (from the second step on) and `start_step()`: move time and the random number streams on
2	people, then every module	`record_states()`: the state at the start of the step
3	custom modules, then demographics	`step()`: custom modules; births add agents, and deaths are requested
4	diseases	`step_state()`: state updates of a disease written as code
5	connectors, networks, interventions	`step()`: cross-module code, then edges, then deliveries
6	diseases	`step()`: transmission, and every transition of a declared disease
7	people	`step_die()`: requested deaths take effect, and every module’s `on_death()` is called
8	people, then every module	`record_events()`: what happened in the step
9	analyzers	`step()`, reading the finished step
10	every module, people, sim	`finish_step()`: dead agents are removed

This is the same order as in Starsim 3. New agents exist before diseases see them, so a child born during a step can be infected in it. The progression of a disease written as code comes before the networks, so that a change of state can change who is eligible for contact (a declared disease takes all its transitions at row 6, so that a delivery sees the same people under every solver). Networks update before interventions, so an intervention sees the final network, and interventions run before transmission, so someone vaccinated today is protected today. Within a kind, modules run in the order given. Three rules follow:

Each module sees what the modules before it did. Infection, testing, and treatment can happen in one step. HIV and syphilis can each change the other’s transmission: the order decides who goes first, and its effect shrinks with the step.
Order inside a declared disease matters little. Transitions out of one state compete in one draw: the total hazard decides who leaves, and each one’s fraction of the total decides where they go, so declaring recovery before death favors neither on average. The few places where declaration order does decide are below.
Nobody is removed mid-step. A death is a request that takes effect at step_die(), so someone who dies during a step can still transmit in it, and arrays never change length partway through.

The plan can be printed, changed, and profiled (Reading the loop). In a compartmental run with nothing that needs a step of its own, the whole run is one call to the kernels (Engines and transmission); sim.describe() says why a run is stepped instead, and require_kernels=True makes that an error.

Per-module timesteps

A module can have its own dt. The loop’s plan is every (time, function) pair of every module on its own timeline, sorted by time and then by the order above, so on the days yearly demographics are due they run at their usual place in a daily sim. Rates and durations are converted with the dt of the module that uses them. describe() shows the schedule; for a daily disease with ss.Births(dt=ss.years(1)) and ss.Deaths(dt=ss.years(1)):

schedule     every day: sir
             every 365 days: births, deaths

How the state is updated within a step

This section is the reference for what happens inside one step; other pages link here. Four questions have separate answers.

1. Competing exits. Exits out of one state that are driven by rates compete: the total rate decides who leaves, and each exit’s share of it decides where they go. In the ODE the flows simply add. The stochastic compartmental solvers draw how many leave a compartment and then split them by share, taking the flows in order of the names of their states and transitions, so the order they are declared in doesn’t change the draws. As agents, each person takes one draw, and the order of declaration sets which part of it goes to which exit: with the same seed, reordering can change which agents take each exit, but not how many do on average. Declaration order does decide a few things as agents: an exact tie between timed exits (two fixed delays of the same length) goes to the one declared first, as does someone reached by two triggers in the same step, and a loop of states is broken at the transition declared first. Exits with p= are decided when people arrive in the state.

2. When the rates are read. The compartmental solvers work out every rate from the state at the start of each substep: the ODE’s own substeps, the chain binomial’s substeps (or the whole step, with ss.solvers.binomial(substeps=1)), tau-leaping’s leaps, or each event under Gillespie. As agents, a declared disease takes its step in two parts:

first the transitions that are not transmission (progression, recovery, death), whose rates all read the counts at the start of the step, so these are simultaneous;
then transmission, which reads who is infectious after those moves, so someone who recovers during a step doesn’t transmit in it.

To first order, the second part is balanced: someone who arrives in a state that transmits has only half a step of its exit rates at their first draw, so each stay counts in transmission for as long as it lasts. Likewise, people who arrive in a state during the first part get half a step of its rates in that step. Both half steps are those of the default, ss.solvers.abm(step='midstep'); with step='wholestep', each rate acts only over whole steps, and with step='events', agents move in continuous time within the step.

3. Cascades. As agents, moves with no wait happen in the step that caused them. After the transitions are taken, anyone who has reached a state whose exit is already due (a trigger whose event has just happened, or a wait that rounds to no steps) moves on, in an order worked out from the graph, until nobody moves (Engines and transmission describes this sweep). The step’s infections are dated at its end and then cascade the same way, so with a zero wait from E to I, someone infected in a step is in I at the next time point. A loop back to the same state waits at least a step. In a compartmental run, a branch with dur=0 goes straight to its destination, and other arrivals leave from the next substep.

4. Between modules. Modules run one after another, in the order of a step, and each sees what the modules before it did, under every solver. A declared disease takes all its transitions in its step() (row 6), so connectors, networks, and interventions see it as it was at the start of the step, and it sees what they changed. A move made from code (an intervention’s move, ss.move(), disease.infect()) is a round of its own, taken when it is made, with its own random draws. A disease written as code does what its code says, changing its states in place in step_state() (row 4) and step() (row 6), as in Starsim 3.

A small example, as agents, of an SEIR with a zero wait from E to I, in one step: the people due to recover leave I (read from the start-of-step counts); transmission then reads who is in I, and infects some people in S; those infections are dated at the end of the step, move to E, and, with no wait, on to I, so they are counted in I at the next time point; then an analyzer (row 9) sees the finished step.

Why it matters: accuracy and the size of the step

Every scheme that takes whole steps is an approximation, with an error that shrinks as the step shrinks. Updating states one after another within a disease, without the half-step balance above, would be the worst choice: someone infected and then recovered in the same step would never transmit, which lowers the effective R0 by a fraction of order γ·dt. This is why the agent engine reads its rates at the start of the step and balances what transmission sees. Between modules, the error is of order dt, like the rest of the step. sim.check_dt() measures every effect of order dt at once; see Solvers for how big it is under each solver, and what to do about it.

Errors, warnings, and notes

The policy is: prefer guessing a default and printing it to raising an error. An error is for when there is no sensible guess: a name that doesn’t exist, lengths that don’t match, a request the solver can’t meet. Every message is meant to be read once, by an epidemiologist or by an AI writing a model for one, so it uses the user’s own names, says what exists instead, and says what to do:

S, I, R = ss.states('S I R')
sir = ss.Disease(infection=S >> I @ 0.05 * I, recovery=I >> R @ ss.days(7), init_prev=0.01, name='sir')
sim = ss.Sim(sir, pars=dict(betta=0.2), verbose=0)
try:
    sim.init()
except KeyError as e:
    print(f'{type(e).__name__}: {e}')

KeyNotFoundError: pars= has a value for "betta", which matches no parameter. Parameters are sir.init_prev, sir.beta, sir.dur_I; did you mean sir.beta?

There are three levels:

Level	Type	What happens
Error	a built-in exception: a `KeyError` for a name that is not found, a `TypeError` for a value of the wrong kind, a `ValueError` for the rest	raises; caught by type, e.g. `except KeyError`
Warning	a `RuntimeWarning`	issued through Python’s warnings, so filters (and `ss.options.warnings='error'`) apply; the run continues
Note	`ss.note()`	a line about something assumed or approximated; kept in `sim.notes` whatever the verbosity, and printed as one line when `verbose` is 1 or more

A note is said once: the same note isn’t kept twice, and notes made while a sim is being initialized are held and printed after the line that starts the run, so they read as part of it. Since the notes are kept, they survive into a saved sim, and a test can check for them.

The few errors and warnings that code catches or filters by kind (ss.CapabilityError, ss.AlreadyRunError, ss.errors.AtBound, and others) are classes of their own, in starsim/internals/errors.py, listed in Profiling and debugging. The tests check each error, warning, and note by its message (tests/devtests/test_errors.py), so a change to a message is a deliberate one.

Advanced: importers and exporters

Every importer reads another framework’s model into one intermediate object, ss.converters.Source: its states, parameters, flows as rate expressions, stratifications, outcomes, and what it couldn’t translate. One function, ss.converters.write(), writes a Source as Starsim code (a declared ss.Disease, with its parameters, states, and transitions as variables, ending with sim = ss.Sim(...)), by the same rules for every format, and ss.converters.build() runs that code to make the sim. So there is one path, from a foreign model to code to a sim, and never from a sim back to code; each framework needs only a reader, and there are no framework-to-framework converters. The code is kept as sim.imported.code, for the user to keep and change. A value the writer can’t write as code, such as a network, is given by the importer as ss.converters.Code('...'), the Python that makes it.

flowchart LR
  F["odin, SBML, EMULSION, SimInf,<br>flepiMoP, summer2, Atomica, camdl,<br>EpiModel, epipack, epydemix, epymorph, MEmilio"] --> R[reader] --> S["ss.converters.Source"] --> W["ss.converters.write()"] --> C["Starsim code"] --> Sim[sim]
  Sim --> E["SBML, odin"]

The hardest part of importing is reading rates. Most frameworks write a flow as a total, such as beta * S * I / N; the importer turns it into a per-person rate on a transition, and says how it read each one. A rate divided by the population total is exposure through the routes (beta * I); one without the division is density-dependent and is kept as a count (beta * I.n); anything else is kept literally, as an expression of counts, which is exact under every solver.

An importer keeps the source’s names, marks what it guessed (as [guessed] in describe()), keeps going when something doesn’t map onto a built-in piece, and lists what it couldn’t translate (ss.errors.ImportIncomplete, also in sim.imported.incomplete). Exporters write SBML and odin from the built model’s compartmental form, and list what they left out or changed (ss.errors.LossyExport); a network, a module written as code, and a function are always left out. A sim made in Starsim is never written back out as Starsim code: its code is the script that made it.

The code is in starsim/converters/: core.py holds Source, write(), and build(), there is one file per framework, and export.py has the exporters. The files for epipack, epydemix, epymorph, MEmilio, and summer2 also run as scripts in those frameworks’ own environments, writing the model as JSON for the importer (EpiModel’s equivalent is dump_epimodel.R). Each importer is tested against its framework’s own examples (tests/devtests/io/). For the modeler’s view, with worked examples, see Import and export.

Limits

A disease written as code is not data. It runs in the same loop, but describe() shows only its class, arrays, and parameters, and it can’t be exported. See Declared modules and modules written as code.
Order between modules is a modeling choice that Starsim won’t make for you. By convention, the module whose accuracy matters most runs last, since it sees the most up-to-date state.