Skip to content

Commit

Permalink
allow to have multiple predictions in one sim.dir
Browse files Browse the repository at this point in the history
  • Loading branch information
@hanase
hanase committed Sep 27, 2024
1 parent 882c754 commit 63c877b
Show file tree
Hide file tree
Showing 5 changed files with 77 additions and 55 deletions.
1 change: 1 addition & 0 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ export(
get.thinned.tfr.mcmc,
get.tfr.prediction,
has.tfr.prediction,
available.tfr.predictions,
get.tfr.convergence.all,
get.tfr.convergence,
get.tfr3.convergence.all,
Expand Down
23 changes: 18 additions & 5 deletions R/get_outputs.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -222,14 +222,23 @@ create.thinned.tfr.mcmc <- function(mcmc.set, thin=1, burnin=0, output.dir=NULL,
return(thin.mcmc.set)
}

has.tfr.prediction <- function(mcmc=NULL, sim.dir=NULL) {
has.tfr.prediction <- function(mcmc=NULL, sim.dir=NULL, subdir = "predictions") {
if (!is.null(mcmc)) sim.dir <- if(is.character(mcmc)) mcmc else mcmc$meta$output.dir
if (is.null(sim.dir)) stop('Either mcmc or directory must be given.')
if(file.exists(file.path(sim.dir, 'predictions', 'prediction.rda'))) return(TRUE)
if(file.exists(file.path(sim.dir, subdir, 'prediction.rda'))) return(TRUE)
return(FALSE)
}

get.tfr.prediction <- function(mcmc=NULL, sim.dir=NULL, mcmc.dir=NULL) {
available.tfr.predictions <- function(mcmc=NULL, sim.dir=NULL, full.names = FALSE){
if (!is.null(mcmc)) sim.dir <- if(is.character(mcmc)) mcmc else mcmc$meta$output.dir
if (is.null(sim.dir)) stop('Either mcmc or directory must be given.')
all.dirs <- list.dirs(sim.dir, recursive = FALSE, full.names = TRUE)
pred.dirs <- names(unlist(sapply(all.dirs, function(x) list.files(x, pattern = "^prediction.rda$"))))
if(length(pred.dirs) == 0 || full.names == TRUE) return(pred.dirs)
return(sapply(pred.dirs, basename, USE.NAMES = FALSE))
}

get.tfr.prediction <- function(mcmc=NULL, sim.dir=NULL, mcmc.dir=NULL, subdir = "predictions") {
############
# Returns an object of class bayesTFR.prediction
# Set mcmc.dir to NA, if the prediction object should not have a pointer
Expand All @@ -238,11 +247,15 @@ get.tfr.prediction <- function(mcmc=NULL, sim.dir=NULL, mcmc.dir=NULL) {
if (!is.null(mcmc))
sim.dir <- if(is.character(mcmc)) mcmc else mcmc$meta$output.dir
if (is.null(sim.dir)) stop('Either mcmc or directory must be given.')
output.dir <- file.path(sim.dir, 'predictions')
output.dir <- file.path(sim.dir, subdir)
pred.file <- file.path(output.dir, 'prediction.rda')
if(!file.exists(pred.file)) {
warning('File ', pred.file, ' does not exist.')
return(NULL)
if(length((alt.preds <- available.tfr.predictions(sim.dir = sim.dir))) > 0){
output.dir <- file.path(sim.dir, alt.preds[1])
pred.file <- file.path(output.dir, 'prediction.rda')
warning('Extracting predictions from ', alt.preds[1])
} else return(NULL)
}
load(file=pred.file)
bayesTFR.prediction$output.directory <- normalizePath(output.dir)
Expand Down
53 changes: 28 additions & 25 deletions R/predict_tfr.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ tfr.predict <- function(mcmc.set=NULL, end.year=2100,
use.tfr3=TRUE, burnin3=2000,
mu=2.1, rho=0.8859, sigmaAR1=0.1016,
min.tfr=0.5, use.correlation=FALSE,
save.as.ascii=0, output.dir = NULL,
save.as.ascii=0, output.dir = NULL, subdir = "predictions",
low.memory=TRUE,
seed=NULL, verbose=TRUE, uncertainty=FALSE, ...) {
if(!is.null(mcmc.set)) {
Expand Down Expand Up @@ -55,7 +55,7 @@ tfr.predict <- function(mcmc.set=NULL, end.year=2100,
start.year=start.year, nr.traj=nr.traj, burnin=burnin, thin=thin, use.tfr3=has.phase3, burnin3=burnin3,
mu=mu, rho=rho, sigmaAR1 = sigmaAR1, min.tfr=min.tfr, use.correlation=use.correlation,
save.as.ascii=save.as.ascii,
output.dir=output.dir, verbose=verbose, uncertainty=uncertainty, ...))
output.dir=output.dir, subdir = subdir, verbose=verbose, uncertainty=uncertainty, ...))
}

.find.burnin.nr.traj.from.diag <- function(diag.list, verbose = FALSE) {
Expand Down Expand Up @@ -88,17 +88,18 @@ get.burnin.nrtraj.from.diagnostics <- function(sim.dir, ...) {
}

tfr.predict.extra <- function(sim.dir=file.path(getwd(), 'bayesTFR.output'),
prediction.dir=sim.dir, countries = NULL, save.as.ascii=0, verbose=TRUE, uncertainty=FALSE,
prediction.dir=sim.dir, subdir = "predictions", countries = NULL, save.as.ascii=0,
verbose=TRUE, uncertainty=FALSE,
all.countries.required = TRUE, use.correlation = NULL) {
# Run prediction for given countries/regions (as codes). If they are not given it will be set to countries
# for which there are MCMC results but no prediction.
# It is to be used after running run.tfr.mcmc.extra
mcmc.set <- get.tfr.mcmc(sim.dir)
if(is.null(mcmc.set))
stop('Error in "sim.dir" argument.')
pred <- get.tfr.prediction(sim.dir=prediction.dir)
pred <- get.tfr.prediction(sim.dir=prediction.dir, subdir = subdir)
if(is.null(pred))
stop('Error in "prediction.dir" argument.')
stop('Error in "sim.dir", "prediction.dir" or/and "subdir" argument. Use available.tfr.predictions() to check on valid predictions directories.')
if(length(setdiff(pred$mcmc.set$meta$regions$country_code, mcmc.set$meta$regions$country_code)) > 0)
stop('Prediction is inconsistent with the mcmc results. Use tfr.predict.')
if(is.null(countries)) {
Expand All @@ -125,7 +126,7 @@ tfr.predict.extra <- function(sim.dir=file.path(getwd(), 'bayesTFR.output'),
use.correlation=use.correlation,
mu=pred$mu, rho=pred$rho, sigmaAR1=pred$sigmaAR1,
min.tfr=pred$min.tfr, countries=countries.idx, save.as.ascii=0, output.dir=prediction.dir,
force.creating.thinned.mcmc=TRUE,
subdir = subdir, force.creating.thinned.mcmc=TRUE,
write.summary.files=FALSE, write.trajectories=TRUE, verbose=verbose,
uncertainty=uncertainty, all.countries.required = all.countries.required)

Expand Down Expand Up @@ -173,7 +174,8 @@ make.tfr.prediction <- function(mcmc.set, start.year=NULL, end.year=2100, replac
use.correlation=FALSE, countries = NULL,
adj.factor1=NA, adj.factor2=0, forceAR1=FALSE,
boost.first.period.in.phase2=TRUE,
save.as.ascii=0, output.dir = NULL, write.summary.files=TRUE,
save.as.ascii=0, output.dir = NULL, subdir = "predictions",
write.summary.files=TRUE,
is.mcmc.set.thinned=FALSE, force.creating.thinned.mcmc=FALSE,
write.trajectories=TRUE,
verbose=verbose, uncertainty=FALSE, all.countries.required = TRUE){
Expand Down Expand Up @@ -213,10 +215,10 @@ make.tfr.prediction <- function(mcmc.set, start.year=NULL, end.year=2100, replac

#setup output directory
if (is.null(output.dir)) output.dir <- meta$output.dir
outdir <- file.path(output.dir, 'predictions')
outdir <- file.path(output.dir, subdir)

if(is.null(countries)) {
if(!replace.output && has.tfr.prediction(sim.dir=output.dir))
if(!replace.output && has.tfr.prediction(sim.dir=output.dir, subdir = subdir))
stop('Prediction in ', outdir,
' already exists.\nSet replace.output=TRUE if you want to overwrite existing projections.')
unlink(outdir, recursive=TRUE)
Expand Down Expand Up @@ -1005,11 +1007,11 @@ do.convert.trajectories <- function(pred, n, output.dir, countries=NULL, verbose


convert.tfr.trajectories <- function(dir=file.path(getwd(), 'bayesTFR.output'),
n=1000, output.dir=NULL,
n=1000, subdir = "predictions", output.dir=NULL,
verbose=FALSE) {
# Converts all trajectory rda files into UN ascii, selecting n trajectories by equal spacing.
if(n <= 0) return()
pred <- get.tfr.prediction(sim.dir=dir)
pred <- get.tfr.prediction(sim.dir=dir, subdir = subdir)
if (is.null(output.dir)) output.dir <- pred$output.directory
if(!file.exists(output.dir)) dir.create(output.dir, recursive=TRUE)
cat('Converting trajectories from', dir, '\n')
Expand All @@ -1028,11 +1030,11 @@ convert.tfr.trajectories <- function(dir=file.path(getwd(), 'bayesTFR.output'),
}

write.projection.summary <- function(dir=file.path(getwd(), 'bayesTFR.output'),
output.dir=NULL, revision=NULL, adjusted=FALSE,
subdir = "predictions", output.dir=NULL, revision=NULL, adjusted=FALSE,
est.uncertainty = FALSE, ...) {
# Writes three prediction summary files, one in a user-friendly format, one in a UN-format,
# and one parameter file.
pred <- get.tfr.prediction(sim.dir=dir)
pred <- get.tfr.prediction(sim.dir=dir, subdir = subdir)
if (is.null(output.dir)) output.dir <- pred$output.directory
if(!file.exists(output.dir)) dir.create(output.dir, recursive=TRUE)
tfr.write.projection.summary.and.parameters(pred, output.dir, revision=revision, adjusted=adjusted,
Expand Down Expand Up @@ -1302,19 +1304,19 @@ get.tfr.shift <- function(country.code, pred) {
return(pred)
}

tfr.median.reset <- function(sim.dir, countries = NULL) {
tfr.median.reset <- function(sim.dir, countries = NULL, ...) {
if(is.null(countries)) {
pred <- get.tfr.prediction(sim.dir)
pred <- get.tfr.prediction(sim.dir, ...)
pred$median.shift <- NULL
store.bayesTFR.prediction(pred)
cat('\nMedians for all countries reset.\n')
} else
for(country in countries) pred <- tfr.median.shift(sim.dir, country, reset=TRUE)
for(country in countries) pred <- tfr.median.shift(sim.dir, country, reset=TRUE, ...)
invisible(pred)
}

tfr.median.shift <- function(sim.dir, country, reset=FALSE, shift=0, from=NULL, to=NULL) {
pred <- get.tfr.prediction(sim.dir)
tfr.median.shift <- function(sim.dir, country, reset=FALSE, shift=0, from=NULL, to=NULL, ...) {
pred <- get.tfr.prediction(sim.dir, ...)
pred <- .bdem.median.shift(pred, type='tfr', country=country, reset=reset,
shift=shift, from=from, to=to)
store.bayesTFR.prediction(pred)
Expand Down Expand Up @@ -1354,16 +1356,17 @@ tfr.median.shift <- function(sim.dir, country, reset=FALSE, shift=0, from=NULL,
return(pred)
}

tfr.median.set <- function(sim.dir, country, values, years=NULL) {
pred <- get.tfr.prediction(sim.dir)
tfr.median.set <- function(sim.dir, country, values, years=NULL, ...) {
pred <- get.tfr.prediction(sim.dir, ...)
pred <- .bdem.median.set(pred, 'tfr', country=country, values=values, years=years)
store.bayesTFR.prediction(pred)
invisible(pred)
}

tfr.median.adjust <- function(sim.dir, countries, factor1=2/3, factor2=1/3, forceAR1=FALSE) {
pred <- get.tfr.prediction(sim.dir)
if (is.null(pred)) stop('No valid prediction in ', sim.dir)
tfr.median.adjust <- function(sim.dir, countries, factor1=2/3, factor2=1/3, forceAR1=FALSE, subdir = "predictions") {
pred <- get.tfr.prediction(sim.dir, subdir = subdir)
if (is.null(pred)) stop('Prediction not found in ', file.path(sim.dir, subdir),
'. Check available.tfr.predictions() and use the subdir argument to set non-standard prediction subdirectory.')
mcmc.set <- pred$mcmc.set
if(is.null(countries)) {
cat('\nNo countries given. Nothing to be done.\n')
Expand All @@ -1382,15 +1385,15 @@ tfr.median.adjust <- function(sim.dir, countries, factor1=2/3, factor2=1/3, forc
use.tfr3=pred$use.tfr3, mcmc3.set=m3.set, burnin3=pred$burnin3,
mu=pred$mu, rho=pred$rho, sigmaAR1=pred$sigmaAR1, min.tfr=pred$min.tfr,
countries=countries.idx, adj.factor1=factor1, adj.factor2=factor2,
forceAR1=forceAR1, save.as.ascii=0, output.dir=sim.dir,
forceAR1=forceAR1, save.as.ascii=0, output.dir=sim.dir, subdir = subdir,
write.summary.files=FALSE, is.mcmc.set.thinned=TRUE,
write.trajectories=FALSE, verbose=FALSE)
new.means <- new.pred$traj.mean.sd[,1,2:dim(new.pred$traj.mean.sd)[3]]
for(icountry in 1:length(countries)) {
tfr.median.set(sim.dir, countries[icountry], new.means[get.country.object(countries[icountry], mcmc.set$meta)$index,])
}
# reload adjusted prediction
invisible(get.tfr.prediction(sim.dir))
invisible(get.tfr.prediction(sim.dir, subdir = subdir))
}

tfr.correlation <- function(meta, cor.pred=NULL, low.coeffs=c(0.11, 0.26, 0.05, 0.09),
Expand Down
4 changes: 2 additions & 2 deletions R/predict_tfr_subnat.R
View file
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
tfr.predict.subnat <- function(countries, my.tfr.file, sim.dir=file.path(getwd(), 'bayesTFR.output'),
end.year=2100, start.year=NULL, output.dir = NULL, annual = NULL,
end.year=2100, start.year=NULL, subdir = "predictions", output.dir = NULL, annual = NULL,
nr.traj=NULL, seed = NULL, min.tfr = 0.5, ar.pars = NULL,
save.as.ascii = 0, verbose = TRUE) {
# Run subnational projections, using the Scale AR(1) model applied to a national bayesTFR simulation
# sim.dir is the world-national simulation. Set output.dir to store results somewhere else.

wpred <- get.tfr.prediction(sim.dir) # contains national projections
wpred <- get.tfr.prediction(sim.dir, subdir = subdir) # contains national projections
wdata <- wpred$tfr_matrix_reconstructed
wmeta <- wpred$mcmc.set$meta
if(!is.null(seed)) set.seed(seed)
Expand Down
51 changes: 28 additions & 23 deletions R/tfr_median_set_new.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,9 +30,9 @@ tfr.shift.estimation.to.wpp <- function(sim.dir, ..., verbose = TRUE){
meta$median.shift.estimation[[as.character(cntry)]] <- merged$shift
}
store.bayesTFR.meta.object(meta, meta$output.dir)
has.predictions <- has.tfr.prediction(sim.dir = sim.dir)
if(has.predictions) { # need to save this also in the collapsed chain of the prediction
pred <- get.tfr.prediction(sim.dir = sim.dir)
predictions <- available.tfr.predictions(sim.dir = sim.dir)
for(pred.subdir in predictions) { # need to save this also in the collapsed chain of the prediction
pred <- get.tfr.prediction(sim.dir = sim.dir, subdir = pred.subdir)
pred$mcmc.set$meta$median.shift.estimation <- meta$median.shift.estimation
store.bayesTFR.prediction(pred)
store.bayesTFR.meta.object(pred$mcmc.set$meta, pred$mcmc.set$meta$output.dir)
Expand All @@ -47,12 +47,13 @@ get.tfr.shift.estimation <- function(country.code, meta)
return(meta$median.shift.estimation[[as.character(country.code)]])
}

tfr.median.set.all <- function(sim.dir, country, values, years=NULL, burnin = 0, thin = 1)
tfr.median.set.all <- function(sim.dir, country, values, years=NULL, burnin = 0, thin = 1,
subdir = "predictions")
{
mcmc.set <- get.tfr.mcmc(sim.dir)
meta <- mcmc.set$meta
country.obj <- get.country.object(country, meta=meta)
has.predictions <- has.tfr.prediction(sim.dir = sim.dir)
has.predictions <- has.tfr.prediction(sim.dir = sim.dir, subdir = subdir)
# browser()
bdem.shift.estimation <- do.call('get.tfr.shift.estimation', list(country.obj$code, meta))
estimation.years <- as.numeric(dimnames(meta$tfr_matrix)[[1]])
Expand All @@ -63,7 +64,7 @@ tfr.median.set.all <- function(sim.dir, country, values, years=NULL, burnin = 0,
years.missing <- is.null(years)
if (has.predictions)
{
pred <- get.tfr.prediction(sim.dir = sim.dir)
pred <- get.tfr.prediction(sim.dir = sim.dir, subdir = subdir)
pred.years <- as.numeric(dimnames(pred$quantiles)[[3]])
bdem.shift <- do.call('get.tfr.shift', list(country.obj$code, pred))
nr.proj <- pred$nr.projections + 1
Expand Down Expand Up @@ -130,8 +131,10 @@ tfr.median.reset.estimation <- function(sim.dir, countries = NULL)
{
mcmc.set <- get.tfr.mcmc(sim.dir)
meta <- mcmc.set$meta
has.prediction <- has.tfr.prediction(sim.dir = sim.dir)
if (has.prediction) pred <- get.tfr.prediction(sim.dir = sim.dir)
has.prediction <- has.tfr.prediction(sim.dir = sim.dir, subdir = subdir)
#if (has.prediction) pred <- get.tfr.prediction(sim.dir = sim.dir)
predictions <- available.tfr.predictions(sim.dir = sim.dir)

has.estimation <- (!is.null(mcmc.set$mcmc.list[[1]]$uncertainty) && mcmc.set$mcmc.list[[1]]$uncertainty)
output <- list()
if (has.estimation && !is.null(meta$median.shift.estimation))
Expand All @@ -146,23 +149,25 @@ tfr.median.reset.estimation <- function(sim.dir, countries = NULL)
}
output[['meta']] <- meta
store.bayesTFR.meta.object(meta, meta$output.dir)
if (has.prediction){
for(pred.subdir in predictions) { # need to save this also in the prediction mcmc object
pred <- get.tfr.prediction(sim.dir = sim.dir, subdir = pred.subdir)
pred$mcmc.set$meta$median.shift.estimation <- meta$median.shift.estimation
store.bayesTFR.prediction(pred)
}
}
if (has.prediction && !is.null(pred$median.shift))
{
if(is.null(countries)) pred$median.shift <- NULL # reset all countries
else {
for (country in countries)
{
country.obj <- get.country.object(country, meta=meta)
pred$median.shift[[as.character(country.obj$code)]] <- NULL
}
}
output[['pred']] <- pred
store.bayesTFR.prediction(pred)
for(pred.subdir in predictions) {
pred <- get.tfr.prediction(sim.dir = sim.dir, subdir = pred.subdir)
if(is.null(pred$median.shift)) next
if(is.null(countries)) pred$median.shift <- NULL # reset all countries
else {
for (country in countries)
{
country.obj <- get.country.object(country, meta=meta)
pred$median.shift[[as.character(country.obj$code)]] <- NULL
}
}
output[['pred']] <- pred # this just returns the last prediction object which is not correct if there are more than one prediction directories
store.bayesTFR.prediction(pred)
}
invisible(output)
}
Expand Down Expand Up @@ -206,8 +211,8 @@ tfr.median.reset.estimation <- function(sim.dir, countries = NULL)
}


tfr.shift.prediction.to.wpp <- function(sim.dir, ...){
pred <- get.tfr.prediction(sim.dir)
tfr.shift.prediction.to.wpp <- function(sim.dir, subdir = "predictions", ...){
pred <- get.tfr.prediction(sim.dir, subdir = subdir)
new.pred <- .do.shift.prediction.to.wpp(pred, wpp.dataset = "tfrprojMed", ...)
store.bayesTFR.prediction(new.pred)
invisible(new.pred)
Expand Down

0 comments on commit 63c877b

Please sign in to comment.