R/addstep.R

In vincnardelli/epidsampler: Simulations for Sampling in Epidemics

#' Add step
#'
#' Simulate one step (i.e. one day) of mobility, social interaction and the evolution of the epidemics with the infection and evolution of the disease for each individual in the map.
#'
#' @param map An epidmap object
#' @param m A positive float. % of moving people
#' @param s A positive integer. Step size
#' @param cn A positive integer. Mean of contacts inside a tile
#' @param cp A positive integer. Mean of number of people for each contact
#' @param im A positive integer. Mean of number of people infected
#' @param tE A positive integer. Number of days of E condition
#' @param tI A positive integer. Number of days of I condition
#' @param tA A positive integer. Number of days of A condition
#' @param ir A positive float. Infected rate I/(A+I)
#' @param cfr A positive float.case Fatality rate
#' @param verbose A logical
#'
#' @return an epidmap object
#' @export
#'
#' @details
#' ## Mobility
#' The contagion mechanism is favoured by people mobility. In this simulation, we assumed that in any moment of time a certain percentage **m** of the population can move between the squares. In this way it is possible to distinguish different epidemic phases such as free-to-move period and lockdown. The commuting during the lockdown period is not only limited by the number of people who move, but also by the extent of their movements. This is a further simulation parameter which is generated by a uniform distribution ranging from **-s** to **s**.
#' ## Social interaction
#'Given the mobility pattern described above, contagion is determined by the social interaction and the contact opportunities. The number of contacts in each square of the grid is assumed to be determined by a random number drawn from a Poisson distribution with parameter, say **cn**, while the number of people involved in the movements is also a Poisson number characterised but a different parameter **cp**. Given these assumptions, a contagion occurs in the following way. If in a meeting it is present at least one asymptomatic or an exposed person, **im** susceptibles will be infected moving in the status of the Exposed.
#' ## Evolution of the epidemics
#' First of all, in order to simulate an artificial population describing the time evolution of an epidemics, we considered a popular model constituted by a system of six differential equations which, in each moment of time, describe six categories of individuals, namely: the susceptibles (S), those exposed to the virus (E), the infected with symptoms (I), those without symptoms (A) and those that are removed from population either because healed (R) or dead (D). This modelling framework is due to the seminal contribution of Hamer (1906), Kermack and McKendrick (1927) and Soper (1929) and it is often referred to as the “SIR model” from the initials of the categories considered in the first simplified formulation: Susceptibles, Infected and Removed.
#'For the data random generation, we assumed that, if infected, a susceptible element of the population (S) will remain in the exposed state (E) for the time **tE**. After that period the subject can become either infected with symptoms (I) with probability **ir** or without (asymptomatic; symbol A) with probability **1-ir**. The asymptomatic will remain infected (and so still able to transmit the virus) for **tA** days. After this period all the asymptomatic will be considered healed and will pass to the category removed (R). In contrast, the infected people showing symptoms will be healed with probability **1-cfr** or die (D) with probability **cfr** (case death rate).
#'
#'
#' @examples
#' \dontrun{
#'set.seed(12345)
#'map <- genmap(n=5, p=c(10, 50))
#'map <- map %>%
#'  addstep(m=0.03, s=2,
#'           cn=4, cp=2,
#'           im=3, tE=5, tA=14,
#'           tI=14, ir=1, cfr=0.15)
#' }
#' @importFrom stats rpois
addstep <- function(map=map, m=NULL, s=NULL,
                    cn=NULL, cp=NULL,
                    im=NULL,
                    tE=NULL, tI=NULL, tA=NULL, ir=NULL,
                    cfr=NULL, verbose=FALSE){
  n <- map$par$n
  if(!is.null(map$par$m) & is.null(m)) m <- map$par$m
  if(!is.null(map$par$s) & is.null(s)) s <- map$par$s
  if(!is.null(map$par$cn) & is.null(cn)) cn <- map$par$cn
  if(!is.null(map$par$cp) & is.null(cp)) cp <- map$par$cp
  if(!is.null(map$par$im) & is.null(im)) im <- map$par$im
  if(!is.null(map$par$tE) & is.null(tE)) tE <- map$par$tE
  if(!is.null(map$par$tI) & is.null(tI)) tI <- map$par$tI
  if(!is.null(map$par$tA) & is.null(tA)) tA <- map$par$tA
  if(!is.null(map$par$ir) & is.null(ir)) ir <- map$par$ir
  if(!is.null(map$par$cfr) & is.null(cfr)) cfr <- map$par$cfr

  t <-map$par$t + 1
  map$par$t <- t



  move <- function(map, m, s){
    lastmap <- map$data[map$data$t == t-1, ]
    can_move <- lastmap$id[lastmap$condition %in% c("S", "E", "A")]
    n_can_move <- round(length(can_move)*m, 0)


    move_p <- sample(can_move, n_can_move)

    move_x <- round(runif(n_can_move, min=-s, max=s), 0)
    move_y <- round(runif(n_can_move, min=-s, max=s), 0)

    check_border <- function(list){
      list[list < 1] <- 1 - list[list < 1]
      list[list > n] <- n - (list[list > n] - n)
      list[list < 1] <- 1
      list[list > n] <- n

      return(list)
    }

    new_x <- check_border(lastmap[move_p, ]$x + move_x)
    new_y <- check_border(lastmap[move_p, ]$y + move_y)

    #t <- map$par$t+1
    lastmap$t <- t
    lastmap[move_p,]$x <- new_x
    lastmap[move_p,]$y <- new_y

    return(lastmap)
  }

  meet <- function(t, cell, cp){
    #contact in one cell

    #cell to xy
    x <- map$map[cell,]$x
    y <- map$map[cell, ]$y


    cell_data <- map$data[map$data$x == x & map$data$y ==y  & map$data$t == t & map$data$condition %in% c("S", "E", "A", "R"),]
    n_people_contact <- rpois(1, cp)

    if(nrow(cell_data) > n_people_contact & nrow(cell_data) > 1){
      sampled <- sample(cell_data$id, n_people_contact)

      contact <- expand.grid(t=t,
                             p1=sampled,
                             p2=sampled)
      contact$c1 <- sapply(contact$p1, function(k) cell_data[cell_data$id == k, ]$condition)
      contact$c2 <- sapply(contact$p2, function(k) cell_data[cell_data$id == k, ]$condition)

      contact <- contact[contact$p1!=contact$p2,]

      if(nrow(contact) > 0){

        # contagion
        contact$new <- FALSE

        condition <- sapply(sampled, function(k) cell_data[cell_data$id == k, ]$condition)
        if(sum(c("E", "A") %in% condition)!=0 & length(sampled[condition == "S"])>0){
          if(length(sampled[condition == "S"])==1){
            new_exposed <- sampled[condition == "S"]
          }else{
            new_exposed <- sample(sampled[condition == "S"], im)
          }

          contact[contact$p1 %in% new_exposed, ]$new <- TRUE

        }
      }


      return(contact)
    }else{
      return(NULL)
    }

  }

  check_stationary_status <- function(vector, t){
    if(length(vector)< t) return(FALSE)
    vector <- vector[(length(vector)-t+1):length(vector)]
    if(verbose) print(vector)
    last_element <- vector[length(vector)]
    return(mean(vector == rep(last_element, length(vector))) == 1)
  }


  cat("simulating t", t, "\n")
  # number of meeting for each cell
  # TODO cn depends on population
  cns <- rep(1:nrow(map$map), rpois(nrow(map$map), cn))
  # structure
  # move
  map$data <- rbind(map$data, move(map, 0.05, 3))
  # update

  # change status
  # if E
  ids <- map$data$id[map$data$t == max(map$data$t) & map$data$condition == "E"]

  for(idx in ids){
    if (check_stationary_status(map$data$condition[map$data$id == idx], tE)){
      map$data$condition[map$data$t == max(map$data$t) & map$data$id == idx] <- sample(c("I", "A"), 1, prob=c(ir, 1-ir))
      if(verbose) print(map$data$condition[map$data$t == max(map$data$t) & map$data$id == idx])
    }
  }

  # if A

  ids <- map$data$id[map$data$t == max(map$data$t) & map$data$condition == "A"]

  for(idx in ids){
    if (check_stationary_status(map$data$condition[map$data$id == idx], tA)){
      map$data$condition[map$data$t == max(map$data$t) & map$data$id == idx] <- "R"
      if(verbose) print(map$data$condition[map$data$t == max(map$data$t) & map$data$id == idx])
    }
  }

  # if I
  ids <- map$data$id[map$data$t == max(map$data$t) & map$data$condition == "I"]

  for(idx in ids){
    if (check_stationary_status(map$data$condition[map$data$id == idx], tI)){
      map$data$condition[map$data$t == max(map$data$t) & map$data$id == idx] <- sample(c("D", "R"), 1, prob=c(cfr, 1-cfr))
      if(verbose) print(map$data$condition[map$data$t == max(map$data$t) & map$data$id == idx])
    }
  }
  # finish I
  # meet
  possible_meet <- purrr::possibly(meet, otherwise = NULL)
  list <- purrr::map(cns, possible_meet, t=t, cp=cp)

  # future::plan(future::multisession)
  # list <- furrr::future_map(cns, possible_meet, t=t, cp=cp)


  meetdf <- purrr::reduce(list, rbind)
  map$contacts <- rbind(map$contacts, meetdf)
  new_exposed <- unique(meetdf[meetdf$new,]$p1)
  if(length(new_exposed) > 0) map$data[map$data$t == t & map$data$id %in% new_exposed,]$condition <- "E"

  map$par$m <- m
  map$par$s <- s
  map$par$t <- t
  map$par$cn <- cn
  map$par$cp <- cp
  map$par$im <- im
  map$par$tE <- tE
  map$par$tI <- tI
  map$par$tA <- tA
  map$par$ir <- ir
  map$par$cfr <- cfr
  return(map)
}