Solutions

These are the solutions to the problems from each of the tutorials. (Note: Some solutions are still pending, but email us if you’re impatient!)

An interactive version of this notebook is available on Google Colab or Binder.

Let’s start with the simplest version of a Starsim model. We’ll make a version of a classic SIR model. Here’s how our code would look:

T1 Solutions

Question 1

Q: To simulate a susceptible-infectious-susceptible (SIS) model instead of SIR, what would we change in the example above?

A: We would simply change 'sir' to 'sis':

import starsim as ss
import sciris as sc

# Define the parameters
pars = sc.objdict( # We use objdict to allow "." access
    n_agents = 10_000,
    networks = sc.objdict(
        type = 'random',
        n_contacts = 10,
    ),
    diseases = sc.objdict(
        type = 'sis', # <-- change this
        init_prev = 0.01,
        beta = 0.05,
    )
)

# Make the sim, run and plot
sim = ss.Sim(pars)
sim.run()
sim.plot()
sim.diseases.sis.plot() # <-- change this

Initializing sim with 10000 agents
  Running 2000.01.01 ( 0/51) (0.00 s)  ———————————————————— 2%
  Running 2010.01.01 (10/51) (0.15 s)  ••••———————————————— 22%
  Running 2020.01.01 (20/51) (0.18 s)  ••••••••———————————— 41%
  Running 2030.01.01 (30/51) (0.21 s)  ••••••••••••———————— 61%
  Running 2040.01.01 (40/51) (0.24 s)  ••••••••••••••••———— 80%
  Running 2050.01.01 (50/51) (0.28 s)  •••••••••••••••••••• 100%

Figure(768x576)

Figure(672x480)

Question 2

Q: How do the results change if we increase/decrease beta?

Increasing beta makes the curves steeper:

pars.diseases.type = 'sir' # Switch back to SIR
pars2 = sc.dcp(pars) # copy to new dictionary
pars2.diseases.beta = 0.10
sim2 = ss.Sim(pars2).run()
sim2.diseases.sir.plot()

Initializing sim with 10000 agents
  Running 2000.01.01 ( 0/51) (0.00 s)  ———————————————————— 2%
  Running 2010.01.01 (10/51) (0.04 s)  ••••———————————————— 22%
  Running 2020.01.01 (20/51) (0.08 s)  ••••••••———————————— 41%
  Running 2030.01.01 (30/51) (0.16 s)  ••••••••••••———————— 61%
  Running 2040.01.01 (40/51) (0.21 s)  ••••••••••••••••———— 80%
  Running 2050.01.01 (50/51) (0.24 s)  •••••••••••••••••••• 100%

Figure(672x480)

Decreasing beta makes the curves shallower:

pars3 = sc.dcp(pars)
pars3.diseases.beta = 0.02
sim3 = ss.Sim(pars3).run()
sim3.diseases.sir.plot()

Initializing sim with 10000 agents
  Running 2000.01.01 ( 0/51) (0.00 s)  ———————————————————— 2%
  Running 2010.01.01 (10/51) (0.09 s)  ••••———————————————— 22%
  Running 2020.01.01 (20/51) (0.12 s)  ••••••••———————————— 41%
  Running 2030.01.01 (30/51) (0.16 s)  ••••••••••••———————— 61%
  Running 2040.01.01 (40/51) (0.19 s)  ••••••••••••••••———— 80%
  Running 2050.01.01 (50/51) (0.22 s)  •••••••••••••••••••• 100%

Figure(672x480)

Question 3

Q: How do the results change if we reduce the number of agents to 200?

We get a similar result as before, except less smooth, since random effects are more important with small numbers of agents:

pars4 = sc.dcp(pars)
pars4.n_agents = 200
sim4 = ss.Sim(pars4).run()
sim4.diseases.sir.plot()

Initializing sim with 200 agents
  Running 2000.01.01 ( 0/51) (0.00 s)  ———————————————————— 2%
  Running 2010.01.01 (10/51) (0.01 s)  ••••———————————————— 22%
  Running 2020.01.01 (20/51) (0.02 s)  ••••••••———————————— 41%
  Running 2030.01.01 (30/51) (0.04 s)  ••••••••••••———————— 61%
  Running 2040.01.01 (40/51) (0.05 s)  ••••••••••••••••———— 80%
  Running 2050.01.01 (50/51) (0.06 s)  •••••••••••••••••••• 100%

Figure(672x480)

T2 Solutions

Question 1

Q: How would you model an outbreak of an SIR-like disease within a refugee camp of 2,000 people? Suppose you were interested in the cumulative number of people who got infected over 1 year - how would you find this out?

The answer obviously depends on the disease parameters. However, we can make some simple assumptions and use cum_infections to determine the total number of infections:

import starsim as ss
import sciris as sc

pars = sc.objdict(
    n_agents = 2_000,
    start = '2025-01-01',
    dur = 365,
    dt = 'day',
    verbose = 1/30, # Print every month
)
sir = ss.SIR(
    dur_inf = ss.days(14),
    beta = ss.perday(0.02),
    init_prev = 0.001,
)
net = ss.RandomNet(n_contacts=4)

sim = ss.Sim(pars, diseases=sir, networks=net)
sim.init()
sim.run()
sim.plot()

answer = sim.results.sir.cum_infections[-1]
print(f'Cumulative infections over one year: {answer}')

Initializing sim with 2000 agents
  Running 2025.01.01 ( 0/366) (0.00 s)  ———————————————————— 0%
  Running 2025.01.31 (30/366) (0.04 s)  •——————————————————— 8%
  Running 2025.03.02 (60/366) (0.08 s)  •••————————————————— 17%
  Running 2025.04.01 (90/366) (0.12 s)  ••••———————————————— 25%
  Running 2025.05.01 (120/366) (0.16 s)  ••••••—————————————— 33%
  Running 2025.05.31 (150/366) (0.20 s)  ••••••••———————————— 41%
  Running 2025.06.30 (180/366) (0.25 s)  •••••••••——————————— 49%
  Running 2025.07.30 (210/366) (0.29 s)  •••••••••••————————— 58%
  Running 2025.08.29 (240/366) (0.33 s)  •••••••••••••——————— 66%
  Running 2025.09.28 (270/366) (0.37 s)  ••••••••••••••—————— 74%
  Running 2025.10.28 (300/366) (0.41 s)  ••••••••••••••••———— 82%
  Running 2025.11.27 (330/366) (0.44 s)  ••••••••••••••••••—— 90%
  Running 2025.12.27 (360/366) (0.48 s)  •••••••••••••••••••— 99%
Figure(768x576)

Cumulative infections over one year: 411.0

Question 2

Solution pending.

T3 Solutions

Question 1

Q: In Niger, the crude birth rate is 45 and the crude death rate is 9. Assuming these rates stay constant, and starting with a total population of 24 million in 2020, how many people will there be in 2040? (You do not need to include any diseases in your model.)

A: We can build our simple demographic model with these parameters, then run it and plot the results:

import starsim as ss
import sciris as sc

pars = sc.objdict(
    start = 2020,
    stop = 2040,
    total_pop = 24e6,
    birth_rate = 45,
    death_rate = 9,
)
sim = ss.Sim(pars)
sim.run()
sim.plot('n_alive')

answer = sim.results.n_alive[-1]/1e6
print(f'Population size in year {pars.stop}: {answer} million')

Initializing sim with 10000 agents
  Running 2020.01.01 ( 0/21) (0.00 s)  ———————————————————— 5%
  Running 2030.01.01 (10/21) (0.01 s)  ••••••••••—————————— 52%
  Running 2040.01.01 (20/21) (0.02 s)  •••••••••••••••••••• 100%

Figure(768x576)

Population size in year 2040: 49.596 million

T4 Solutions

Question 1

Solution pending.

T5 Solutions

Question 1

Solution pending.

T6 Solutions

Question 1

Q: If we change the disease from SIR to SIS and set coverage to 100%, what minimum efficacy of vaccine is required to eradicate the disease by 2050?

A: There are many ways we could solve this, including with formal numerical optimization packages. However, since we are only varying a single parameter, we can also just use a simple binay search or grid search. This solution illustrates both approaches.

import numpy as np
import sciris as sc
import starsim as ss

pars = dict(
    n_agents = 5_000,
    birth_rate = 20,
    death_rate = 15,
    networks = dict(
        type = 'random',
        n_contacts = 4
    ),
    diseases = dict(
        type = 'sis',
        dur_inf = 10,
        beta = 0.1,
    ),
    verbose = False,
)

class sis_vaccine(ss.Vx):
    """ A simple vaccine against "SIS" """
    def __init__(self, efficacy=1.0, **kwargs):
        super().__init__()
        self.define_pars(efficacy=efficacy)
        self.update_pars(**kwargs)
        return

    def administer(self, people, uids):
        people.sis.rel_sus[uids] *= 1-self.pars.efficacy
        return
    
def run_sim(efficacy):
    """ Run a simulation with a given vaccine efficacy """
    # Create the vaccine product
    product = sis_vaccine(efficacy=efficacy)

    # Create the intervention
    intervention = ss.routine_vx(
        start_year=2015, # Begin vaccination in 2015
        prob=1.0,        # 100% coverage
        product=product  # Use the SIS vaccine
    )

    # Now create two sims: a baseline sim and one with the intervention
    sim = ss.Sim(pars=pars, interventions=intervention)
    sim.run()
    return sim

def objective(efficacy, penalty=10, boolean=False, verbose=False):
    """ Calculate the objective from the simulation """
    sim = run_sim(efficacy=efficacy)
    transmission = sim.results.sis.new_infections[-1] > 0
    if boolean:
        return not transmission
    else:
        loss = efficacy + penalty*transmission
        if verbose:
            print(f'Trial: {efficacy=}, {transmission=}, {loss=}')
        return loss

def grid_search(n=5, reps=2):
    """ Perform a grid search over the objective function """
    sc.heading('Performing grid search ...')
    lb = 0 # Lower bound for efficacy
    ub = 1 # Upper bound for efficacy
    for rep in range(reps):
        print(f'Grid search {rep+1} of {reps}...')
        efficacy = np.linspace(lb, ub, n)
        transmission = sc.parallelize(objective, efficacy, boolean=True)
        lb = efficacy[sc.findlast(transmission, False)]
        ub = efficacy[sc.findfirst(transmission, True)]
        print(f'  Trials: {dict(zip(efficacy, transmission))}')
        print(f'  Results: lower={lb}, upper={ub}')
    mid = (lb+ub)/2
    print(sc.ansi.bold(f'Result: {mid}'))
    return mid, lb, ub

def auto_search(efficacy=1.0):
    """ Perform automatic search """
    sc.heading('Performing automatic search...')
    out = sc.asd(objective, x=efficacy, xmin=0, xmax=1, maxiters=10, verbose=True)
    print(sc.ansi.bold(f'Result: {out.x}'))
    return out

# Run both optimizations
mid, lb, ub = grid_search()
out = auto_search()

——————————————————————————

Performing grid search ...

——————————————————————————



Grid search 1 of 2...

  Trials: {0.0: False, 0.25: False, 0.5: True, 0.75: True, 1.0: True}

  Results: lower=0.25, upper=0.5

Grid search 2 of 2...

  Trials: {0.25: False, 0.3125: False, 0.375: False, 0.4375: False, 0.5: True}

  Results: lower=0.4375, upper=0.5

Result: 0.46875





——————————————————————————————

Performing automatic search...

——————————————————————————————



     step 1 (0.2 s) ++ (orig:array([1.]) | best:array([1.]) | new:array([0.9]) | diff:array([-0.1]))

     step 2 (0.4 s) ++ (orig:array([1.]) | best:array([0.9]) | new:array([0.7]) | diff:array([-0.2]))

     step 3 (0.6 s) -- (orig:array([1.]) | best:array([0.7]) | new:array([10.3]) | diff:array([9.6]))

     step 4 (0.8 s) -- (orig:array([1.]) | best:array([0.7]) | new:array([0.8]) | diff:array([0.1]))

     step 5 (1.0 s) ++ (orig:array([1.]) | best:array([0.7]) | new:array([0.5]) | diff:array([-0.2]))

     step 6 (1.2 s) -- (orig:array([1.]) | best:array([0.5]) | new:array([10.1]) | diff:array([9.6]))

     step 7 (1.5 s) -- (orig:array([1.]) | best:array([0.5]) | new:array([0.55]) | diff:array([0.05]))

     step 8 (1.7 s) -- (orig:array([1.]) | best:array([0.5]) | new:array([10.3]) | diff:array([9.8]))

     step 9 (1.8 s) -- (orig:array([1.]) | best:array([0.5]) | new:array([10.4]) | diff:array([9.9]))

     step 10 (2.0 s) -- (orig:array([1.]) | best:array([0.5]) | new:array([10.45]) | diff:array([9.95]))

===  Maximum iterations reached (10 steps, orig: 1.000 | best: 0.5000 | ratio: 0.49999999999999994) ===

Result: 0.49999999999999994