#' branch
#' Simulates a continuous-time time-inhomogenous markov branching process using the specified parameters. Uses C++ code for faster simulation.
#'
#' @param model the \code{process_model} object representing the process being simulates
#' @param params the vector of parameters for which we are simulating the model
#' @param time_obs the vector of times at which to record the process state
#' @param reps the number of replicates to simulate
#' @param silent if true, verbose output will be shown. Default: false
#' @param keep if true, the temporary comma-separated file generated while simulating while be kept. if false, it will be deleted. Default: false
#' @param seed seed for the random number generator. If NULL, will use computer clock to set a random seed
#'
#' @export
branch <- function(model, params, init_pop, time_obs, reps, silent = FALSE, keep = FALSE, seed = NULL){
if(class(model) != "estipop_process_model"){
stop("model must be a process_model object!")
}
if((!is.numeric(params) && !is.null(params)) || !is.numeric(init_pop) || !is.numeric(time_obs) || !is.numeric(reps)){
stop("all time, population, and parameter inputs must be numeric!")
}
if(length(init_pop) != model$ntypes){
stop("init_pop and model must have same number of types ")
}
if(any(init_pop < 0) || reps <= 0){
stop("population must be nonnegative and reps must be positive!")
}
if(any(time_obs < 0)){
stop("all observation times must be nonnegative.")
}
if(!is.logical(silent) || !is.logical(keep)){
stop("parameters slient and keep should be logical!")
}
#modify the transition list with metadata about whether the rate is constant
timedep <- F
for(i in 1:length(model$transition_list)){
if(is_const(model$transition_list[[i]]$rate$exp)){
#evaluate constant rates
model$transition_list[[i]]$rate <- eval(model$transition_list[[i]]$rate$exp, list(params = params))
model$transition_list[[i]]$type <- 1
}
else{
timedep <- T
fname <- paste("custom_rate_", digest::digest(deparse(model$transition_list[[i]]$rate),"md5"), sep="")
create_timedep_template(model$transition_list[[i]]$rate$exp, params, paste(fname, ".cpp", sep=""))
compile_timedep( paste(fname, ".cpp", sep=""))
model$transition_list[[i]]$rate <- c(paste(fname, ".so", sep=""), "rate")
model$transition_list[[i]]$type <- 2
}
}
f <- R.utils::getAbsolutePath(tempfile(pattern = paste("system_", format(Sys.time(), "%d-%m-%Y-%H%M%S"), "_", sep = ""), fileext = ".csv", tmpdir = getwd()))
if(timedep){
if(is.null(seed)){
timeDepBranch(time_obs, reps, f, init_pop, model$transition_list, stops = NULL, silent)
} else {
timeDepBranch(time_obs, reps, f, init_pop, model$transition_list, stops = NULL, silent, seed)
}
} else {
if(is.null(seed)){
gmbp3(time_obs, reps, f, init_pop, model$transition_list, stops = NULL, silent)
} else {
gmbp3(time_obs, reps, f, init_pop, model$transition_list, stops = NULL, silent, seed)
}
}
res <- read.csv(f, header = F)
names(res)[1:2] <- c("rep","time")
if(!keep){
file.remove(f)
}
#remove the temporary dynamic libraries
for(trans in model$transition_list){
if(trans$type == 2){
fname = .pop(trans$rate[1],".so")
file.remove(trans$rate[1])
file.remove(paste(fname, ".cpp", sep=""))
file.remove(paste(fname, ".h", sep=""))
file.remove(paste(fname, ".o", sep=""))
file.remove(paste(fname, ".cpp.backup", sep=""))
}
}
res <- data.frame(res)
names(res) <- c("rep","time",paste("type", 1:model$ntypes, sep=""))
return(res)
}
#' branch_approx
#' Approximately simulated a time-inhomogenous Markov branching process by sampling from the asymptotic multivariate normal distribution.
#'
#' @param model the \code{process_model} object representing the process being simulates
#' @param params the vector of parameters for which we are simulating the model
#' @param time_obs the vector of times at which to record the process state
#' @param reps the number of replicates to simulate
#'
#' @export
branch_approx = function(model, params, init_pop, time_obs, reps){
if(class(model) != "estipop_process_model"){
stop("model must be a process_model object!")
}
if((!is.numeric(params) && !is.null(params)) || !is.numeric(init_pop) || !is.numeric(time_obs) || !is.numeric(reps)){
stop("all time, population, and parameter inputs must be numeric!")
}
if(length(init_pop) != model$ntypes){
stop("init_pop and model must have same number of types ")
}
if(any(init_pop < 0) || reps <= 0){
stop("population must be nonnegative and reps must be positive!")
}
if(any(time_obs < 0)){
stop("all observation times must be nonnegative.")
}
if(class(model) != "estipop_process_model"){
stop("model must be a process_model object!")
}
if(length(init_pop) != model$ntypes){
stop("init_pop and model must have same number of types ")
}
if(!is.numeric(params) || !is.numeric(init_pop) || !is.numeric(time_obs) || !is.numeric(reps)){
stop("all time, population, and parameter inputs must be numeric!")
}
ntype = model$ntypes
time_obs <- sort(unique(c(0, time_obs)))
res <- matrix(rep(0, reps*(ntype+2)*length(time_obs)), ncol = ntype + 2)
mom <- moments(model, params, time_obs[1:(length(time_obs)-1)], time_obs[-1])
for(r in 1:reps){
curr_pop <- init_pop
res[(r-1)*length(time_obs)+ 1,] <- c(time = 0, rep = r, curr_pop)
for(k in 2:length(time_obs)){
desol <- as.matrix(mom[(mom$tf == time_obs[k]), -c(1,2)])
m_mat <- matrix(desol[1:ntype**2],nrow= ntype)
dt_mat <- matrix(desol[(ntype**2 + 1):(ntype**3 + ntype**2)],nrow = ntype) #second moment array
mu_vec <- t(m_mat)%*%curr_pop #final mean population vector
sigma_mat <- matrix(curr_pop%*%dt_mat,c(ntype,ntype*ntype)) #final population covariance matrix
for(i in 1:ntype){
for(j in 1:ntype){
sigma_mat[i,j] <- sigma_mat[i,j] - curr_pop%*%(m_mat[,i]*m_mat[,j])
}
}
curr_pop <- round(MASS:::mvrnorm(n = 1, mu_vec, sigma_mat))
res[(r-1)*length(time_obs) + k,] <- c(time = time_obs[k], rep = r, curr_pop)
}
res <- data.frame(res)
names(res) <- c("time","rep",paste("type", 1:ntype, sep=""))
}
return(res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.