How to test - SIR script example

sir_script.R

# parameters
N <- 1000      # network size
avk <- 7       # desired average degree
gamma <- 0.25  # recovery rate
tau <- 0.27    # per-edge infection rate
max_time <- 50 # maximum simulation length
i0 <- 10       # initial number of infected individuals

# build erdos-renyi random-graph
network <- matrix(0L, nrow = N, ncol = N)
p <- avk / (N - 1)
for (i in 1:(N - 1)) {
  for (j in (i + 1):N) {
    if (stats::runif(1) < p){
      network[i, j] <- 1L
      network[j, i] <- 1L
    }
  }
}

# choose some nodes to infect at random
initial_infected <- sample(1:N, i0, replace = FALSE)

# set the initial status of nodes
status <- integer(N)
status[initial_infected] <- 1

# calculate the initial rates for each node
rate <- integer(N)
rate[initial_infected] <- gamma
for (i in initial_infected) {
  susceptible_contacts <- which((network[, i] == 1) & (status == 0))
  rate[susceptible_contacts] <- rate[susceptible_contacts] + tau
}

# preallocate for output (R thing)
times <- rep(NA_real_, 2 * N)
INF <- rep(NA_integer_, 2 * N)
SUS <- rep(NA_integer_, 2 * N)

# prep for simulation
t <- 0   # time
i <- 1   # index
INF[i] <- i0
SUS[i] <- N - i0
times[i] <- 0

# begin simulation
while((t < max_time) && (INF[i] > 0)) {

  total_rate <- sum(rate)
  cumulative_rate <- cumsum(rate)

  # calculate what node is changing status
  event <- match(TRUE, cumulative_rate > (runif(1) * total_rate))

  # calculate timestep to the event and update time
  timestep <- -log(runif(1)) / total_rate
  t <- t + timestep
  times[i + 1] <- t

  # update status of nodes, S -> I (0 -> 1) I -> R (1 -> 2)
  status[event] <- status[event] + 1

  # update rates
  susceptible_contacts <- which((network[, event] == 1) & (status == 0))
  if (status[event] == 1) {
    # S -> I so can recover with rate gamma
    rate[event] <- gamma

    # record infected and susceptible
    INF[i + 1] <- INF[i] + 1
    SUS[i + 1] <- SUS[i] - 1

    # susceptible contacts gain an I neighbour
    rate[susceptible_contacts] <- rate[susceptible_contacts] + tau
  } else {
    # I -> R so no infectious pressure
    rate[event] <- 0

    # susceptible contacts lose an I neighbour
    rate[susceptible_contacts] <- rate[susceptible_contacts] - tau

    # record infected and susceptible
    INF[i + 1] <- INF[i] - 1
    SUS[i + 1] <- SUS[i]
  }

  # increment index
  i <- i + 1
}