Exercise: Python vs R

OverviewPythonR

Here we have provided SIR ODE model implementations in Python and in R. Each script runs several scenarios and produces a plot of infection prevalence for each scenario.

You can download each script to debug on your computer:

• sir_ode.py
• sir_ode.R

#!/usr/bin/env python3

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp


def sir_rhs(time, state, popn, beta, gamma):
    """
    The right-hand side for the vanilla SIR compartmental model.
    """
    s_to_i = beta * state[1] * state[0] / popn  # beta * I(t) * S(t) / N
    i_to_r = gamma * state[1]  # gamma * I(t)
    return [-s_to_i, s_to_i - i_to_r, i_to_r]


def run_model(settings):
    """
    Return the SIR ODE solution for the given model settings.
    """
    # Define the time span and evaluation times.
    sim_days = settings['sim_days']
    time_span = [0, sim_days]
    times = np.linspace(0, sim_days, num=sim_days + 1)
    # Define the initial state.
    popn = settings['population']
    exposures = settings['exposures']
    initial_state = [popn - exposures, exposures, 0]
    # Define the SIR parameters.
    params = (popn, settings['beta'], settings['gamma'])
    # Return the daily number of people in S, I, and R.
    return solve_ivp(
        sir_rhs, time_span, initial_state, t_eval=times, args=params
    )


def default_settings():
    """
    The default model settings.
    """
    return {
        'sim_days': 20,
        'population': 100,
        'exposures': 2,
        'beta': 1.0,
        'gamma': 0.5,
    }


def run_model_scaled_beta(settings, scale):
    """
    Adjust the value of ``beta`` before running the model.
    """
    settings['beta'] = scale * settings['beta']
    return run_model(settings)


def run_model_scaled_gamma(settings, scale):
    """
    Adjust the value of ``gamma`` before running the model.
    """
    settings['gamma'] = scale * settings['gamma']
    return run_model(settings)


def plot_prevalence_time_series(solutions):
    """
    Plot daily prevalence of infectious individuals for one or more scenarios.
    """
    fig, axes = plt.subplots(
        constrained_layout=True,
        nrows=len(solutions),
        sharex=True,
        sharey=True,
    )
    for ix, (scenario_name, solution) in enumerate(solutions.items()):
        ax = axes[ix]
        ax.title.set_text(scenario_name)
        ax.plot(solution.y[1], label='I')
        ax.set_xticks([0, 5, 10, 15, 20])
    # Save the figure.
    png_file = 'sir_ode_python.png'
    fig.savefig(png_file, format='png', metadata={'Software': None})


def demonstration():
    settings = default_settings()
    default_scenario = run_model(settings)
    scaled_beta_scenario = run_model_scaled_beta(settings, scale=1.5)
    scaled_gamma_scenario = run_model_scaled_gamma(settings, scale=0.7)

    plot_prevalence_time_series(
        {
            'Default': default_scenario,
            'Scaled Beta': scaled_beta_scenario,
            'Scaled Gamma': scaled_gamma_scenario,
        }
    )


if __name__ == '__main__':
    demonstration()

#!/usr/bin/env -S Rscript --vanilla

library(deSolve)
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(ggplot2))


# The right-hand side for the vanilla SIR compartmental model.
sir_rhs <- function(time, state, params) {
  s_to_i <- params$beta * state["I"] * state["S"] / params$popn
  i_to_r <- params$gamma * state["I"]
  list(c(-s_to_i, s_to_i - i_to_r, i_to_r))
}


# Return the SIR ODE solution for the given model settings.
run_model <- function(settings) {
  # Define the time span and evaluation times.
  times <- seq(0, settings$sim_days)
  # Define the initial state.
  popn <- settings$population
  exposures <- settings$exposures
  initial_state <- c(S = popn - exposures, I = exposures, R = 0)
  # Define the SIR parameters.
  params <- list(
    popn = settings$population,
    beta = settings$beta,
    gamma = settings$gamma
  )
  # Return the daily number of people in S, I, and R.
  ode(initial_state, times, sir_rhs, params)
}


# The default model settings.
default_settings <- function() {
  list(
    sim_days = 20,
    population = 100,
    exposures = 2,
    beta = 1.0,
    gamma = 0.5
  )
}


# Adjust the value of ``beta`` before running the model.
run_model_scaled_beta <- function(settings, scale) {
  settings$beta <- scale  * settings$beta
  run_model(settings)
}


# Adjust the value of ``gamma`` before running the model.
run_model_scaled_gamma <- function(settings, scale) {
  settings$gamma <- scale * settings$gamma
  run_model(settings)
}


# Plot daily prevalence of infectious individuals for one or more scenarios.
plot_prevalence_time_series <- function(solutions) {
  df <- lapply(
    names(solutions),
    function(name) {
      solutions[[name]] |>
        as.data.frame() |>
        mutate(scenario = name)
    }
  ) |>
    bind_rows() |>
    mutate(
      scenario = factor(scenario, levels = names(solutions), ordered = TRUE)
    )
  fig <- ggplot() +
    geom_line(aes(time, I), df) +
    xlab(NULL) +
    scale_y_continuous(
      name = NULL,
      limits = c(0, 40),
      breaks = c(0, 20, 40)
    ) +
    facet_wrap(~ scenario, ncol = 1) +
    theme_bw(base_size = 10) +
    theme(
      strip.background = element_blank(),
      panel.grid = element_blank(),
    )
  png_file <- "sir_ode_r.png"
  ggsave(png_file, fig, width = 640, height = 480, units = "px", dpi = 150)
}


demonstration <- function() {
  settings <- default_settings()
  default_scenario <- run_model(settings)
  scaled_beta_scenario <- run_model_scaled_beta(settings, scale=1.5)
  scaled_gamma_scenario <- run_model_scaled_gamma(settings, scale=0.7)

  plot_prevalence_time_series(
    list(
      Default = default_scenario,
      `Scaled Beta` = scaled_beta_scenario,
      `Scaled Gamma` = scaled_gamma_scenario
    )
  )
}

demonstration()

The model outputs differ!

Here are prevalence time-series plots produced by each script:

Python plotR plot

Model outputs for the Python script.

Model outputs for the R script.

Interactive debugger sessions

If your editor supports running a debugger, use this feature! See these examples for RStudio, PyCharm, Spyder, and VS Code.

Some initial thoughts ...

• Is it obvious whether one of the figures is correct and the other is wrong?

• The sir_rhs() functions in the two scripts appear to be equivalent — but are they?

• The default_settings() functions appear to be equivalent — but are they?

• The run_model_scaled_beta() and run_model_scaled_gamma() functions also appear to be equivalent.

• Where might you begin looking?