R/CR_model.R

Defines functions CR_model

Documented in CR_model

#' Simulate consumer-resource dynamics
#' @description 
#' This function takes a named vector as initial community, simulates based
#' on CR_ode function, and reports the output result matrix. 
#' @name CR_model
#' @param community the initial state of species and resource in a community. A named vector.
#' @param time_limit time limit for ode. 
#' @param time_step time step for ode.
#' @return a deSolve matrix 
#' @examples
#' # Load example parameters
#' data(example_parameter)
#' 
#' # Build regional species pool
#' pool_build()
#' 
#' # Create a community, a named vector
#' set.seed(1)
#' species <- sample(1:S_pool, S, replace = FALSE) %>% sort()
#' community <- 
#'   c(setNames(c(rep(1, P)), paste0("R", sprintf("%03d", 1:P))),
#'     setNames(c(rep(.1, length(species))), paste0("X", sprintf("%05d", species))))
#' 
#' # Run ode
#' result <- CR_model(community = community, time_limit = 100, time_step = 1)
#' 
#' # Plot
#' result %>%
#'  as.tibble() %>%
#'  gather("index", "value", -1) %>%
#'  ggplot(aes(x = time, y = value, color = index, lty = substr(index, 1, 1))) +
#'  geom_line() + 
#'  theme(legend.position="none")
#' @export

# Main function for running consumer-resource model
CR_model <- function (
  community = community, # Named vector with biomass as elements and consumer or resource in names.
  time_limit = 100,  # Time limit
  time_step = 1
) {
  
  # Consumer 
  sp <- grep("X", names(community), value = T) %>% substr(2, 6) %>% as.numeric()
  
  # Resource
  re <- grep("R", names(community), value = T) %>% substr(2, 4) %>% as.numeric()
  
  # Time sequence 
  time <- seq(0, time_limit, by = time_step)
  
  # Parameters: a named vector
  parameters <- list(sp = sp, re = re, J = J, Rin = Rin, P = P, u = u, w = w, m = m, stoi = stoi)
  
  # Initial condition: a named vector
  state <- community

  # Run ode
  result <- ode(y = state, times = time, func = CR_ode, parms = parameters)
  
  # Melt result
  result <- 
    result %>%
    as.data.frame() %>%
    melt(id.var = "time") %>% 
    as.tibble()
  
  #
  return(result)
}
Chang-Yu-Chang/MigrationCommunity documentation built on April 22, 2019, 7:36 p.m.