Calculating derived outputs

In the previous example, we saw how to create and run a compartmental model.

This example shows you how you can request more detailed outputs from the model, in addition to just the compartment sizes. Summer supports the calculation of “derived outputs”: these are additional outputs that are calculated from either:

• the model compartment sizes for each timestep; or

• the model flow rates at each timestep

There are several different types of derived outputs that will be presented in this example:

• `Flow output </api/model.html#summer.model.CompartmentalModel.request_output_for_flow>`__: A sum of one or more flow rates at each timestep

• `Compartment output </api/model.html#summer.model.CompartmentalModel.request_output_for_compartments>`__: A sum of one or more compartment sizes at each timestep

• `Aggregate output </api/model.html#summer.model.CompartmentalModel.request_aggregate_output>`__: An aggregate of other derived outputs

• `Cumulative output </api/model.html#summer.model.CompartmentalModel.request_cumulative_output>`__: A cumulative sum of another derived output

• `Function output </api/model.html#summer.model.CompartmentalModel.request_function_output>`__: A pure function of other derived outputs

To start, let’s define some utility functions to create a SIR model that is similar to the one from the last example:

import numpy as np
import matplotlib.pyplot as plt
from summer import CompartmentalModel

def build_model():
    """Returns a new SIR model"""
    model = CompartmentalModel(
        times=[0, 20],
        compartments=["S", "I", "R"],
        infectious_compartments=["I"],
        timestep=0.1,
    )
    model.set_initial_population(distribution={"S": 990, "I": 10})
    model.add_infection_frequency_flow(name="infection", contact_rate=2, source="S", dest="I")
    model.add_transition_flow(name="recovery", fractional_rate=1/3, source="I", dest="R")
    model.add_death_flow(name="infection_death", death_rate=0.05, source="I")
    return model

def plot_output(model, name, title):
    """Plot a derived output for a model"""
    fig, ax = plt.subplots(1, 1, figsize=(12, 6), dpi=120)
    ax.plot(model.times, model.derived_outputs[name])
    ax.set_title(title)
    ax.set_xlabel("Days")
    ax.set_ylabel(name)
    start, end = ax.get_xlim()
    ax.xaxis.set_ticks(np.arange(start + 1, end, 5))
    plt.show()

# Force NumPy to format printed arrays nicely.
np.set_printoptions(formatter={'all': lambda f: f"{f:0.2f}"})

Let’s quickly visualize what the compartments are doing over time:

model = build_model()
model.run()
fig, ax = plt.subplots(1, 1, figsize=(12, 6), dpi=120)

# Add each compartment to the plot.
for i in range(model.outputs.shape[1]):
    ax.plot(model.times, model.outputs.T[i])

ax.set_title("SIR Model Outputs")
ax.set_xlabel("Days")
ax.set_ylabel("Compartment size")
ax.legend(["S", "I", "R"])
start, end = ax.get_xlim()
ax.xaxis.set_ticks(np.arange(start + 1, end, 5))
plt.show()

Requesting derived outputs

We can ask a model to calculate extra outputs that are derived from the compartment sizes and flow rates. For example, we might want to ask the model to track the number of people who died from infection per timestep.

# Create a model (see above)
model = build_model()

# Request that the model calculate a derived output when it is run.
model.request_output_for_flow(name="deaths", flow_name="infection_death")

Now when we run the model, it will create an NumPy array of infection deaths that we can access via model.derived_outputs.

# Run the model
model.run()

# View the derived outputs dictionary that we calculated when `run()` was called.
print("Derived outputs:\n", model.derived_outputs)

Derived outputs:
 {'deaths': array([0.05, 0.05, 0.06, 0.07, 0.09, 0.10, 0.12, 0.14, 0.16, 0.19, 0.22,
       0.26, 0.30, 0.35, 0.40, 0.46, 0.53, 0.60, 0.69, 0.78, 0.88, 0.99,
       1.10, 1.22, 1.35, 1.47, 1.60, 1.72, 1.84, 1.96, 2.06, 2.16, 2.24,
       2.32, 2.38, 2.43, 2.46, 2.49, 2.50, 2.50, 2.49, 2.48, 2.45, 2.42,
       2.39, 2.35, 2.30, 2.26, 2.21, 2.16, 2.10, 2.05, 2.00, 1.94, 1.88,
       1.83, 1.78, 1.72, 1.67, 1.62, 1.57, 1.51, 1.47, 1.42, 1.37, 1.32,
       1.28, 1.24, 1.20, 1.15, 1.11, 1.08, 1.04, 1.00, 0.97, 0.93, 0.90,
       0.87, 0.84, 0.81, 0.78, 0.75, 0.73, 0.70, 0.67, 0.65, 0.63, 0.60,
       0.58, 0.56, 0.54, 0.52, 0.50, 0.49, 0.47, 0.45, 0.43, 0.42, 0.40,
       0.39, 0.37, 0.36, 0.35, 0.34, 0.32, 0.31, 0.30, 0.29, 0.28, 0.27,
       0.26, 0.25, 0.24, 0.23, 0.22, 0.21, 0.21, 0.20, 0.19, 0.19, 0.18,
       0.17, 0.17, 0.16, 0.15, 0.15, 0.14, 0.14, 0.13, 0.13, 0.12, 0.12,
       0.11, 0.11, 0.11, 0.10, 0.10, 0.09, 0.09, 0.09, 0.08, 0.08, 0.08,
       0.08, 0.07, 0.07, 0.07, 0.07, 0.06, 0.06, 0.06, 0.06, 0.05, 0.05,
       0.05, 0.05, 0.05, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.03,
       0.03, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03, 0.02, 0.02, 0.02,
       0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02,
       0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,
       0.01, 0.01, 0.01])}

Flow outputs

A `flow output </api/model.html#summer.model.CompartmentalModel.request_output_for_flow>`__ tracks a set of requested flow rates for each timestep. These requests can also select flows between particular strata in a stratified model (see later examples).

model = build_model()

# Request that the 'infection_death' flow is tracked as a derived output named 'deaths'.
model.request_output_for_flow(name="deaths", flow_name="infection_death")

model.run()
plot_output(model, "deaths", "Infection deaths per timestep")

Cumulative outputs

You can use a `cumulative output </model.html#summer.model.CompartmentalModel.request_cumulative_output>`__ to request that the model tracks the cumulative sum of other derived outputs over time. For example, let’s track total infection deaths and the total people recovered:

model = build_model()
model.request_output_for_flow(name="deaths", flow_name="infection_death")

# Request that the 'deaths' derived output is accumulated into 'deaths_cumulative'.
model.request_cumulative_output(name="deaths_cumulative", source="deaths")

model.run()
plot_output(model, "deaths_cumulative", "Cumulative infection deaths")

Compartment outputs

A `compartment output </model.html#summer.model.CompartmentalModel.request_output_for_compartments>`__ tracks the sum of one or more compartments at each timestep. These requests can also select compartments for particular strata in a stratified model (see later examples).

model = build_model()

# Request that the S and R compartment sizes are combined into 'uninfected'.
model.request_output_for_compartments(name="uninfected", compartments=["S", "R"])

model.run()
plot_output(model, "uninfected", "Number of currently uninfected people")

Aggregate outputs

You can use an `aggregate output </model.html#summer.model.CompartmentalModel.request_aggregate_output>`__ to request an aggregate of other derived outputs.

model = build_model()

# Track some flows.
model.request_output_for_flow(name="deaths", flow_name="infection_death")
model.request_output_for_flow(name="recoveries", flow_name="recovery")

# Accumulate the flows.
model.request_cumulative_output(name="deaths_cumulative", source="deaths")
model.request_cumulative_output(name="recoveries_cumulative", source="recoveries")

# Aggregate 'deaths_cumulative' and 'recovered_cumulative' into a single output.
model.request_aggregate_output(
    name="dead_or_recovered_cumulative",
    sources=["deaths_cumulative", "recoveries_cumulative"]
)

model.run()
plot_output(model, "dead_or_recovered_cumulative", "Cumulative dead or recovered")
# (In this simple model, this could be also easily be tracked as the complement of the the susceptible population.)

Function outputs

You can use `function outputs </model.html#summer.model.CompartmentalModel.request_function_output>`__ to calculate outputs with a user-defined function. This function takes the requested sources (NumPy arrays) as inputs and should return a NumPy array of the same size.

For example, here we request a calculation that gets us the prevalence of the disease.

model = build_model()

# Track the number of infectious people as a derived output.
# Here we use `save_results=False` because we don't need to use this later.
model.request_output_for_compartments(name="count_infectious", compartments=["I"], save_results=False)

# Track the total population as a derived output.
model.request_output_for_compartments(name="total_population", compartments=["S", "R"], save_results=False)

# Define a function that we will use to calculate prevalence.
# The arguments (infectious, total) will be NumPy arrays.
def get_prevalence(infectious, total):
    return infectious / total

# Request a function output, using `get_prevalence`.
# The sources will map to the function arguments.
model.request_function_output(
    name="prevalence",
    sources=["count_infectious", "total_population"],
    func=get_prevalence,
)

model.run()
plot_output(model, "prevalence", "Prevalence of disease")

Summary

That’s it for now, now you know how to:

• Request derived outputs

• Chain and combine derived outputs

• Access and visualize the derived outputs

A detailed API reference of the CompartmentalModel class can be found here