R/run_model.R

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)

}
kenkellner/jagsUI2 documentation built on July 5, 2019, 9:38 a.m.