Basic Simulator

Now that we are more familiar with the concept of RUMEs, let’s run some simulations.

BasicSimulator is the main tool for the job. It runs one forward simulation at a time and produces one simulation output which records how the model’s state evolved over the duration.

Starting with the same RUME we’ve been using…

import numpy as np
from epymorph.kit import *
from epymorph.adrio import us_tiger

rume = SingleStrataRUME.build(
    ipm=ipm.SIRS(),
    mm=mm.Centroids(),
    init=init.SingleLocation(location=0, seed_size=100),
    scope=StateScope.in_states(["AZ", "CO", "NM", "UT"], year=2020),
    time_frame=TimeFrame.rangex("2020-01-01", "2021-01-01"),
    params={
        "beta": 0.3,
        "gamma": 1/5,
        "xi": 1/90,
        "population": [7_151_502, 5_773_714, 2_117_522, 3_271_616],
        "centroid": us_tiger.InternalPoint(),
    },
)

Running a Simulation

We can construct a simulator passing in this RUME:

sim = BasicSimulator(rume)

And use sim’s run() method to execute the simulation and produce a single simulation output, sometimes called a “realization”:

out = sim.run()

Easy as that!

Overriding Params

The run method gives us a few options worth mentioning. First, we can override some or all of the RUME parameters for the sake of a single realization with the params parameter. For instance, we could quickly explore the effects of different beta values without having to rebuild our RUME.

sim.run(
    params={
        "beta": 0.4,
    },
)

sim.run(
    params={
        "beta": 0.5,
    },
)

This params dictionary follows the same rules as the params dictionary we provided to the RUME, so we don’t lose any flexibility by specifying values this way.

Controlling Randomness

Normally each run is subject to random variations by design. But if we provide a random number generator with a specific seed, we’ll get the same results every time. The rng_factory parameter accepts a function of arity zero (no parameters) that returns a numpy.random.Generator instance.

# default_rng from the epymorph.kit import works
sim.run(rng_factory=default_rng(42))


# or we could define our own function
from numpy.random import Generator, SFC64

def my_rng():
    return Generator(SFC64(42))

sim.run(rng_factory=my_rng)


# or use Python's lambda syntax
sim.run(rng_factory=lambda: Generator(SFC64(42)))

Next we’ll briefly inspect the simulation output that BasicSimulator produces.

a brief introduction to Intra-population Models (IPMs)

A full description of IPMs is included in a later chapter, but it will help to have a basic understanding before continuing. IPMs are epymorph’s way of expressing the transmission and disease mechanics of the simulation. They are compartmental models, where the population is divided into groups representing which disease state its members are in at a given point in time.

The SIRS model has three states:

Susceptible: individuals may become infected if they come into contact with an infectious individual,
Infectious: individuals are currently infected and may spread infection to others, and
Recovered: individuals have recovered from an infection and currently have some immunity against catching the disease again.

When an individual transitions from one state to another, this is an event. Becoming infected, transitioning from susceptible to infectious, is the S-to-I event or \(S \to I\). There are two more events: \(I \to R\) and \(R \to S\). Other transitions are strictly not allowed by the model structure, like going from infectious directly back to susceptible.

All you need to know for now is that these compartments and events exist, and that their ordering is critical to the interpretation of simulation results. When in doubt, we can inspect the ipm object to print this information.

print("compartments:", [str(c.name) for c in rume.ipm.compartments])
print("events:", [str(e.name) for e in rume.ipm.events])

compartments: ['S', 'I', 'R']
events: ['S → I', 'I → R', 'R → S']

Simulation Output

Initial Conditions

Our out object (an instance of the Output class) contains a bunch of neat data. For instance, we can look at the initial conditions of the simulation:

out.initial

array([[7151402,     100,       0],
       [5773714,       0,       0],
       [2117522,       0,       0],
       [3271616,       0,       0]])

The API doc for Output.initial tells us this is an array of shape (N,C) – N for geo nodes and C for IPM compartments. So knowing that, we can now read this array. The first row is for the first location, Arizona, the second for Colorado, and so on. And the columns are the IPM compartments \(S\), \(I\), and \(R\). Just as we would expect, we started with 100 infected individuals in Arizona, and all other individuals are susceptible. No one starts out as recovered – which in epidemiology terms means our population is naive to this hypothetical pathogen.

Compartments

So if that’s how the simulation started, how did things go from there? Imagine that we advance our simulation by one time-step then record this same info. Here we go:

out.compartments[0, :, :]

array([[7151149,     107,       3],
       [5773524,       0,       0],
       [2117855,       0,       0],
       [3271716,       0,       0]])

There’s a few more infected folks and some recovered. Blink forward again:

out.compartments[1, :, :]

array([[7151376,     113,      13],
       [5773714,       0,       0],
       [2117522,       0,       0],
       [3271616,       0,       0]])

Even more folks become infected and some more recover. We get a sense that the simulation is evolving.

The compartments array tracks the time-series data of how many individuals are in each IPM compartment at each simulation location. It’s a (T,N,C)-shaped numpy array, and I’ve used numpy’s indexing features to look at one time-slice at a time. Since this is just a numpy array, we can inspect a few of its critical aspects:

cs = out.compartments
print(f"{type(cs)=}")
print(f"{cs.dtype=}")
print(f"{cs.shape=}")

type(cs)=<class 'numpy.ndarray'>
cs.dtype=dtype('int64')
cs.shape=(732, 4, 3)

And we can use the array to perform interesting calculations, like what was the maximum number of people infected at one time in each US State throughout the simulation? (As we saw, the \(I\) compartment is second in the IPM order, which corresponds to index 1 in Python.)

out.compartments[:, :, 1].max(axis=0)

array([505714, 409055, 149016, 231564])

Events

The compartment values (sometimes called state variables) of our simulation is only one view however. Looking at the above we can’t tell exactly how many people caught the disease during that second time step. The number of infectious individuals changed, but that is a combination of people becoming infected as well as people recovering – there’s an in-flux and out-flux from the \(I\) compartment. In this case we can’t figure it out from looking at the other compartments either, because each of those are also subject to in-flux and out-flux. If you’re just counting the number of people in a room, you get the same number regardless if five leave and five enter or if six leave and six enter.

So we also need to know exactly how many individuals transitioned between each of the compartments at every time step. epymorph records this in the events array.

out.events[1, :, :]

array([[16, 10,  0],
       [ 0,  0,  0],
       [ 0,  0,  0],
       [ 0,  0,  0]])

The events array is (T,N,E)-shaped, where E is the number of possible events in our model. So in this case, we get a row per US State and a column per event: \(S{\to}I\), \(I{\to}R\), and \(R{\to}S\). Now we have the complete picture of how the simulation state changed at every time step.

Data Aggregation Mode

However there’s one more critical detail to keep in mind about how epymorph stores event and compartment data. How should we associate the data to locations? If someone is traveling out of state and they get the flu, does that infection event count towards their home location or the one they’re visiting?

Well, both! Or rather, you might care about either or both of these spatial aggregations depending on the analysis you’re trying to do with the results. epymorph keeps track of both of these datasets for you. So the Output class has events and compartments, as mentioned above, but it also has visit_events, visit_compartments, home_events, and home_compartments.

t = 594  # I picked a time step to demonstrate the difference

# aggregated by visiting location
out.visit_events[t, :, :]

array([[11861, 12035,  8197],
       [11520, 14740,  9268],
       [ 1273,  1189,  2207],
       [ 1498,  1509,  3816]])

# aggregated by home location
out.home_events[t, :, :]

array([[11862, 12039,  8197],
       [11521, 14742,  9271],
       [ 1271,  1184,  2205],
       [ 1498,  1508,  3815]])

And just to make the difference more obvious:

out.visit_events[t, :, :] - out.home_events[t, :, :]

array([[-1, -4,  0],
       [-1, -2, -3],
       [ 2,  5,  2],
       [ 0,  1,  1]])

You can of course access these properties directly. As a convenience, Output objects have a data_mode setting which sets the assumed default data mode. When you call events or compartments, you get whichever arrays correspond to the Output object’s data mode. (Which saves you from having to specify the mode everywhere you use an Output.) The “by visit” mode is the default. If you want a copy of your Output object that uses a specific mode, use the data_by_visit or data_by_mode property.

out_2 = out.data_by_visit
np.array_equal(out_2.events, out.visit_events)

True

out_3 = out.data_by_home
np.array_equal(out_3.events, out.home_events)

True

But regardless of which mode an Output object is in, it still contains both datasets. To add some programming detail, data_by_visit and data_by_home create shallow copies of the Output object; all copies contain references to the same underlying data arrays.

To restate the above in another way, just in case the terminology is a little misleading, “by visit” is not meant to suggest that people are only counted if they are away from their normal home location. Everyone gets counted; it’s just that they count towards where they’re at “right now”, instead of where they came from originally. Similarly, “by home” also counts everyone; it’s just that they count towards where they originated regardless of where they are “right now”.

Other Utilities

Simulation metadata

Aside from the simulation results data, the Output object also holds onto useful metadata about the simulation that produced the output. Most importantly, we can recall our RUME with out.rume and then further query the RUME as needed.

# TimeFrame objects print as ISO-8601 date ranges and duration
print(out.rume.time_frame)

2020-01-01/2020-12-31 (366D)

Output as a DataFrame

If we would prefer to work with our data as a Pandas DataFrame, that’s also supported. We get a table where the time and geo node axes are in “long” format, while the compartments and events are “wide”.

df = out.dataframe

print(f"{type(df)=}")

df.head() # This is a big dataframe, so just display the head of it.

type(df)=<class 'pandas.core.frame.DataFrame'>

	tick	date	node	S	I	R	S → I	I → R
0	0	2020-01-01	04	7151149	107	3	10	3
1	0	2020-01-01	08	5773524	0	0	0	0
2	0	2020-01-01	35	2117855	0	0	0	0
3	0	2020-01-01	49	3271716	0	0	0	0
4	1	2020-01-01	04	7151376	113	13	16	10

Output ticks-to-days conversion

out.ticks_in_days produces an array that converts the simulation time ticks into their equivalent in fractional days. In this simulation, the first simulation step covers one-third of a day while the next covers the remaining two-thirds of a day, and so on repeatedly.

# using np.printoptions() just for nicer formatting...
with np.printoptions(precision=3, floatmode="fixed"):
    print(out.ticks_in_days[0:8])

[0.333 1.000 1.333 2.000 2.333 3.000 3.333 4.000]

Events per day

A common need is to determine how many events happened per day, regardless of how many tau steps are used in the simulation. out.events_per_day is a convenience method to calculate this for you:

# Show the sum of events by location for day-index 123
out.events_per_day[123, :, :]

array([[47801, 59976, 34009],
       [    0,     0,     0],
       [  512,   391,    26],
       [  190,   126,    10]])

Compartments per day

It doesn’t make sense to sum compartment values over time. Compartment counts are a kind of census: the headcount of individuals at a particular point in time. If my fruit bowl contained 5 apples when I looked in the morning and 4 when I checked again in the evening, it would be nonsense to say I had 9 apples.

Instead, if we want to reduce the time axis – going from tau steps to days – for compartments like we did for events, it makes more sense to pick which tau step during the day we care about. For example, I might be interested in the compartment values at the end of each day. In this example, I could use numpy slice operations to yield daily compartments like this:

compartments_daily = out.compartments[1::2, :, :]

compartments_daily[123, :, :]

array([[3780957,  289748, 3080797],
       [5773712,       0,       2],
       [2112387,    1811,    3324],
       [3269857,     683,    1076]])

You can learn more about numpy slicing in their docs. In this case, the slice 1::2 means:

start from time index 1 (skip the first value),
go all the way to the end of the array, and
use a step size of 2 (take every other value).

My simulation uses 2 tau steps, so the second value of every day represents the end-of-day counts (data are recorded after each tau step is applied).