run_model <- function(run_dat){
quiet <- run_dat$quiet
if(quiet|run_dat$n_cores>1) options("jags.pb"="none")
#Initialize output object
out <- list()
#Start modules
load_modules(run_dat$modules)
#Start factories ## NOT IMPLEMENTED YET!
#Get chain ID for printing progress (if necessary)
ch_id <- NULL
if(!is.null(run_dat$this_chain)){
ch_id <- paste0('[chain ',run_dat$this_chain,'] ')
}
#Initialize model object-----------------------------------------------------
if(is.null(run_dat$model_object)){
if(!quiet) cat(paste0(ch_id,'Initializing model\n'))
run_dat$model_object <- init_model(model_file=run_dat$model_file,
data=run_dat$data, inits=run_dat$inits,
n_chains=run_dat$n_chains)
}
#----------------------------------------------------------------------------
#Adaptive phase--------------------------------------------------------------
if(!quiet) cat(paste0(ch_id,'Adapting\n'))
adapt_check <- adapt_model(run_dat$model_object, n_iter=run_dat$n_adapt)
if(any(!adapt_check)){
warning(paste0(ch_id,"JAGS reports adaptation was incomplete"))
}
#----------------------------------------------------------------------------
#Warm-up phase---------------------------------------------------------------
st_time <- Sys.time()
if ( run_dat$n_warmup > 0 ){
if ( run_dat$n_cores == 1 | run_dat$n_warmup < 10 | quiet ){
#Update in one chunk if single-core or few iterations
if(!quiet) cat(paste0(ch_id,'Warm-up\n'))
update_model(run_dat$model_object, n_iter=run_dat$n_warmup)
} else {
#Run in ~10% chunks to show progress
chunks <- get_chunks(run_dat$n_warmup)
for (i in 1:length(chunks)){
update_model(run_dat$model_object, n_iter=chunks[i])
tl <- time_left(st_time, sum(chunks[1:i]), run_dat$n_iter)
pct <- round(sum(chunks[1:i])/run_dat$n_warmup*100)
if(!quiet) cat(paste0(ch_id,'Warm-up ',pct,'%, ',tl,'\n'))
}
}
}
#----------------------------------------------------------------------------
#Sample from posterior-------------------------------------------------------
n_saved <- run_dat$n_iter - run_dat$n_warmup
#Common function for generating samples
get_samples <- function(n_saved){
rjags::coda.samples(model=run_dat$model_object,
variable.names=run_dat$params,
n.iter=n_saved, thin=run_dat$n_thin)
}
#Get samples
if ( run_dat$n_cores == 1 | n_saved < 10 | quiet ){
#Get all samples at once if single-core or few iterations
if (!quiet) cat(paste0(ch_id,'Sampling posterior\n'))
out$samples <- get_samples(n_saved)
} else {
#Run in ~10% chunks to show progress
chunks <- get_chunks(n_saved)
sample_chunks <- vector(length=length(chunks),"list")
for (i in 1:length(chunks)){
sample_chunks[[i]] <- get_samples(chunks[i])
tl <- time_left(st_time, run_dat$n_warmup + sum(chunks[1:i]),
run_dat$n_iter)
pct <- round(sum(chunks[1:i])/n_saved*100)
if(!quiet) cat(paste0(ch_id,'Sampling ', pct,'%, ',tl,'\n'))
}
out$samples <- Reduce(comb_mcmc_list,sample_chunks) #Combine chunks
}
#----------------------------------------------------------------------------
#Save model object
out$model <- run_dat$model_object
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.