inst/scripts/recalc_pdf.R

In JGCRI/gcamland: GCAM Land Allocation Module

#' Recalculate the posterior PDFs for a set of saved scenarios
#'
#' Load a set of scenario-info files generated by an ensemble run,
#' recalculate their posterior probability densities, and optionally save the
#' results.
#'
#' Normally the reason you would do this would be to change the prior
#' and/or likelihood functions used in the analysis, or to add those
#' calculations to scenario outputs where they were never calculated in the
#' first place.  Either way, it only makes sense to do this with uniformly
#' sampled, posterior-weighted samples.  If you do this with Markov Chain Monte
#' Carlo samples, the recorded likelihood and posterior probability values will
#' change, but the sampling will still reflect the old probability density.
#'
#' tldr: \strong{Do not use this function on the results of a Markov Chain Monte
#' Carlo analysis!}  If you do you will get bogus results.
#'
#' The code assumes that the input files have names following
#' \code{\link{run_ensemble}}'s convention:
#' \code{scenario-info-<serialnumber>.rds}.  There is currently no way provided
#' to override this convention.
#'
#' There are several options for saving the output.  \code{TRUE} causes the
#' original files to be overwritten, while \code{NA} (the default) causes the
#' new files to be written into a subdirectory of the original input directory:
#' \code{<inputdir>/revised}.  This subdirectory is created if necessary, but
#' there is no check to see if files in that directory would be overwritten.  If
#' the save option is set to \code{FALSE}, then no output is written.  This
#' option is useful if you just want to get the list of updated ScenarioInfo
#' structures returned by the function.
#'
#' Finally, the results will be calculated in parallel using the
#' \code{doParallel} package.  See the documentation for that package for
#' details on setting up a parallel run.
#'
#' @param inputdir Directory containing the \code{scenario-info*.rds} files.
#' @param saveopt Option for saving (or not) the output.  See Details section.
#' @param lpdf Log-PDF function.
#' @param lprior Log-prior function.
#' @return List of modified ScenarioInfo structures
#' @export
#' @keywords internal
recalc_pdf <- function(inputdir, saveopt=NA, lpdf=get_lpdf(1), lprior=uniform_prior) {

    fn <- NULL                          # silence package checker

    ## Output files to process
    sifiles <- list.files(inputdir, 'scenario-info.*\\.rds')

    ## index start and stop position in filenames
    indx_strt <- nchar('scenario-info-') + 1
    indx_stop <- indx_strt + 5

    if(is.na(saveopt)) {
        outdir <- file.path(inputdir, 'revised')
        if(!dir.exists(outdir)) {
            dir.create(outdir)
        }
    }

    foreach(fn = sifiles) %dopar% {
        file <- file.path(inputdir,fn)
        si <- readRDS(file)

        si_revised <- run_bayes(si, lpdf=lpdf, lprior=lprior)

        if(is.na(saveopt)) {
            outfile <- file.path(outdir, fn)
        }
        else {
            outfile <- file
        }

        if(is.na(saveopt) || saveopt==TRUE) {
            saveRDS(si_revised, outfile)
        }

        si_revised
    }
}