R/simulation_general_functions.R

In tjtnew/simnetR: Epidemic Simulations on Networks

#' Initially Infectious Individuals
#'
#' \code{init_startnodes} returns an integer array of randomly selected values
#' within a given range.
#'
#' This functon is a wrapper around \code{sample} to generate \code{i0} random
#' integers in the interval \code{1:10}.  \code{init_startnodes(N, i0)} is
#' equivalent to \code{sample(1:10, i0, replace = FALSE)}.
#'
#' @param N Number of vertices / nodes in the network.
#' @param i0 Number of initially infectious.
#' @keywords internal
#' @return Integer array of selected individuals.
init_startnodes <- function(N, i0) {
  sample(1:N, i0, replace = FALSE)
}



#' Initial Status of Individuals in the Network
#'
#' \code{init_status} returns a binary array of the initial infectious status
#' of individuals in the network (0 = susceptible, 1 = infectious).
#'
#' @param N Number of vertices / nodes in the network.
#' @param startnodes Initially infectious individuals.
#' @keywords internal
#' @return Integer array of size \code{N} representing binary infectious status.
init_status <- function(N, startnodes) {
  status <- integer(N)
  status[startnodes] <- 1
  return(status)
}



#' Susceptible Contacts of an Individual
#'
#' \code{sus_contacts_matrix} returns a vector of an individuals susceptible
#' neighbours.
#'
#' @param edge_vector The edge vector for individual (taken from adj matrix)
#' @param status The infectious status of all individuals.
#' @keywords internal
#' @return Integer array of susceptible neighbours.
sus_contacts_matrix <- function(edge_vector, status) {
  which((edge_vector == 1) & (status == 0))
}



#' Susceptible Contacts of an Individual
#'
#' \code{sus_contacts_list} returns a vector of an individuals susceptible
#' neighbours.
#'
#' @param contacts The contacts for individual (taken from adj list)
#' @param status The infectious status of all individuals.
#' @keywords internal
#' @return Integer array of susceptible neighbours.
sus_contacts_list <- function(contacts, status) {
  contacts[status[contacts] == 0]
}





#' Infectious Contacts of an Individual
#'
#' \code{inf_contacts_matrix} returns a vector of an individuals infectious
#' neighbours.
#'
#' @param edge_vector The edge vector for individual (taken from adj matrix)
#' @param status The infectious status of all individuals.
#' @keywords internal
#' @return Integer array of infectious neighbours.
inf_contacts_matrix <- function(edge_vector, status) {
  which((edge_vector == 1) & (status == 1))
}


#' Infectious contacts of an Individual
#'
#' \code{inf_contacts_list} returns a vector of an individuals susceptible
#' neighbours.
#'
#' @param contacts The contacts for individual (taken from adj list)
#' @param status The infectious status of all individuals.
#' @keywords internal
#' @return Integer array of susceptible neighbours.
inf_contacts_list <- function(contacts, status) {
  contacts[status[contacts] == 1]
}



#' What Node is Changing Status
#'
#' \code{which_event} returns the node whose status is changing (from S -> I or
#' from I -> S).
#'
#' @param rate The current rate for each nodes
#' @keywords internal
#' @return Integer representing which node the current event is changing
which_event <- function(rate) {
  # total rate
  total <- sum(rate)

  # cumulative rate vector
  cum <- cumsum(rate)

  # find which node the event is happening to (1st non-zero entry)
  event <- match(TRUE, (cum > (stats::runif(1) * total)))

  return(event)
}


#' What is the timestep to the next event
#' \code{step_update} calculates the time to next event
#'
#' @param rate The current rate of nodes in graph
#' @keywords internal
#' @return time to next event
step_update <- function(rate) {
  # draw the timestep based upon the total rate: the higher the rate, the
  # smaller the timestep
  return(-log(stats::runif(1)) / sum(rate))
}