R/get_outputs.R

In PPgp/bayesTFR: Bayesian Fertility Projection

get.tfr.mcmc <- function(sim.dir=file.path(getwd(), 'bayesTFR.output'), chain.ids=NULL, low.memory=TRUE,
							burnin=0, verbose=FALSE) {
	############
	# Returns an object of class bayesTFR.mcmc.set
	############
	mcmc.file.path <- file.path(sim.dir, 'bayesTFR.mcmc.meta.rda')
	if(!file.exists(mcmc.file.path)) {
		warning('File ', mcmc.file.path, ' does not exist.')
		return(NULL)
	}
	load(file=mcmc.file.path)
	bayesTFR.mcmc.meta$output.dir <- normalizePath(sim.dir)
	if(is.null(bayesTFR.mcmc.meta$annual.simulation)) bayesTFR.mcmc.meta$annual.simulation <- FALSE
	if (is.null(chain.ids)) {
		mc.dirs.short <- list.files(sim.dir, pattern='^mc[0-9]+', full.names=FALSE)
		chain.ids <- as.integer(substring(mc.dirs.short, 3))
	} else {
		mc.dirs.short <- paste('mc', chain.ids, sep='')
	}
	ord.idx <- order(chain.ids)
	mc.dirs.short <- mc.dirs.short[ord.idx]
	chain.ids <- chain.ids[ord.idx]
	mcmc.chains <- list()
	counter<-1
	for (imc.d in chain.ids) {
		if (verbose)
			cat('Loading chain', imc.d, 'from disk. ')
		bayesTFR.mcmc <- local({
			load(file=file.path(sim.dir, mc.dirs.short[counter], 'bayesTFR.mcmc.rda'))
			bayesTFR.mcmc})
		mc <- c(bayesTFR.mcmc, list(meta=bayesTFR.mcmc.meta))
		class(mc) <- class(bayesTFR.mcmc)
		if (!low.memory) { # load full mcmc traces
			th.burnin <- get.thinned.burnin(mc, burnin)
			mc$traces <- load.tfr.parameter.traces.all(mc, burnin=th.burnin)
			mc$traces.burnin <- th.burnin
		} else { # traces will be loaded as they are needed
			mc$traces <- 0
			mc$traces.burnin <- 0
		}
		mc$output.dir <- mc.dirs.short[counter]
		if(is.null(mc$uncertainty)) mc$uncertainty <- FALSE
		if (verbose)
			cat('(mcmc.list[[', counter, ']]).\n')
		mcmc.chains[[counter]] <- mc
		counter <- counter+1
	}
	names(mcmc.chains) <- chain.ids
	return(structure(list(meta=bayesTFR.mcmc.meta, 
				mcmc.list=mcmc.chains), class='bayesTFR.mcmc.set'))
}

get.tfr3.mcmc <- function(sim.dir=file.path(getwd(), 'bayesTFR.output'), ...) {
    if(!has.tfr3.mcmc(sim.dir)) stop("Phase III has not been simulated.")
	parent.mc <- get.tfr.mcmc(sim.dir)
	mc <- get.tfr.mcmc(file.path(sim.dir, 'phaseIII'), ...)
	mc$meta$parent <- parent.mc$meta
	mc$meta$regions <- parent.mc$meta$regions
	mc$meta$annual.simulation <- NULL # this needs to be removed as the right value will be taken from meta$parent
	if(length(mc$mcmc.list) > 0) {
		for(chain in 1:length(mc$mcmc.list)) {
			mc$mcmc.list[[chain]]$meta <- mc$meta
		}
	}
	return(mc)					
}

has.tfr.mcmc <- function(sim.dir) {
	return(file.exists(file.path(sim.dir, 'bayesTFR.mcmc.meta.rda')))
}

has.tfr3.mcmc <- function(sim.dir) {
	return(has.tfr.mcmc(file.path(sim.dir, 'phaseIII')))
}

tfr.mcmc <- function(mcmc.set, chain.id) return (mcmc.set$mcmc.list[[chain.id]])

tfr.mcmc.list <- function(mcmc.set, chain.ids=NULL) {
	if(is.null(chain.ids)) return(mcmc.set$mcmc.list)
	return(mcmc.set$mcmc.list[chain.ids])
}

get.thinned.tfr.mcmc <- function(mcmc.set, thin=1, burnin=0) {
	dir.name <- file.path(mcmc.set$meta$output.dir, paste('thinned_mcmc', thin, burnin, sep='_'))
	if(file.exists(dir.name)) return(get.tfr.mcmc(dir.name))
	return(NULL)
}

create.thinned.tfr.mcmc <- function(mcmc.set, thin=1, burnin=0, output.dir=NULL, verbose=TRUE, uncertainty=FALSE,
                                    update.with.countries = NULL) {
	#Return a thinned mcmc.set object with burnin removed and all chanins collapsed into one
    mcthin <- max(sapply(mcmc.set$mcmc.list, function(x) x$thin))
	thin <- max(c(thin, mcthin))
	meta <- mcmc.set$meta
	total.iter <- get.stored.mcmc.length(mcmc.set$mcmc.list, burnin=burnin)
	total.iter <- min(total.iter, get.stored.mcmc.length.extra(mcmc.set$meta, update.with.countries, 
	                                                            nr.chains = length(mcmc.set$mcmc.list), burnin = burnin), na.rm = TRUE)
	meta$is.thinned <- TRUE
	meta$parent.iter <- get.total.iterations(mcmc.set$mcmc.list, burnin)
	meta$parent.meta <- mcmc.set$meta
	meta$nr.chains <- 1
	
	if(verbose) cat('\nStoring thinned mcmc:')
	# store the meta object
	meta$output.dir <- file.path(
		if(is.null(output.dir)) meta$output.dir else output.dir, 
			paste('thinned_mcmc', thin, burnin, sep='_'))
	if(!file.exists(meta$output.dir)) 
		dir.create(meta$output.dir, recursive=TRUE)
	
	if(!is.null(update.with.countries)) meta$one.use.only <- TRUE # thinned traces might get out of sync with other countries; thus recompute next time
	
	store.bayesTFR.meta.object(meta, meta$output.dir)
	
	thin.index <- if(thin > mcthin) unique(round(seq(1, total.iter, by=thin/mcthin))) else 1:total.iter
	nr.points <- length(thin.index)
	
	#create one collapsed mcmc
	thinned.mcmc <- mcmc.set$mcmc.list[[1]]
	thinned.mcmc$meta <- meta
	thinned.mcmc$thin <- 1
	thinned.mcmc$id <- 1
	thinned.mcmc$traces <- 0
	thinned.mcmc$length <- nr.points
	thinned.mcmc$finished.iter <- nr.points
	thinned.mcmc$compression.type <- meta$compression.type
	thinned.mcmc$uncertainty <- uncertainty
	thinned.mcmc$output.dir <- 'mc1'	
	outdir.thin.mcmc <- file.path(meta$output.dir, 'mc1')
	if(!file.exists(outdir.thin.mcmc)) dir.create(outdir.thin.mcmc)

	store.bayesTFR.object(thinned.mcmc, outdir.thin.mcmc)
	
	if(verbose) cat('\nStoring country-independent parameters ...')
	for (par in tfr.parameter.names(trans=FALSE, meta = meta)) {
		values <- get.tfr.parameter.traces(mcmc.set$mcmc.list, par, burnin,
											thinning.index=thin.index)
		write.values.into.file.cindep(par, values, outdir.thin.mcmc, compression.type=thinned.mcmc$compression.type)
	}
	if (!is.null(mcmc.set$meta$ar.phase2) && mcmc.set$meta$ar.phase2) 
	{
	  values <- get.tfr.parameter.traces(mcmc.set$mcmc.list, 'rho_phase2', burnin,
	                                     thinning.index=thin.index)
	  write.values.into.file.cindep('rho_phase2', values, outdir.thin.mcmc, compression.type=thinned.mcmc$compression.type)
	}
	if(verbose) cat('done.\nStoring country-specific parameters ...')
	par.names.cs <- tfr.parameter.names.cs(trans=FALSE)
	if (uncertainty) par.names.cs <- c(par.names.cs, 'tfr')
	cntries <- if(!is.null(update.with.countries)) update.with.countries else 1:mcmc.set$meta$nr_countries
	for (country in cntries){
		country.obj <- get.country.object(country, mcmc.set$meta, index=TRUE)
		if (country.obj$index %in% mcmc.set$meta$extra)
		{
		  nr.points.cs <- floor((mcmc.set$meta$extra_iter[country.obj$index]-burnin) / mcmc.set$meta$extra_thin[country.obj$index])
		  nr.points.cs <- nr.points.cs * length(mcmc.set$mcmc.list)
		  if (nr.points.cs >= nr.points) thin.index.cs <- round(seq(1, nr.points.cs, length.out = nr.points))
		  else 
		  {
		    warning('Country ',country.obj$index,' has less stored iterations. Samples are used repeatedly.\n')
		    thin.index.cs <- rep(1:nr.points.cs, ceiling(nr.points/nr.points.cs))[1:nr.points]
		  }
		}
		for (par in par.names.cs) {
		    if(par != "tfr" && ! country.obj$index %in% mcmc.set$meta$id_DL) next
		  if (country.obj$index %in% mcmc.set$meta$extra)
		  {
		    values <- get.tfr.parameter.traces.cs(mcmc.set$mcmc.list, country.obj, par, 
		                                          burnin=burnin, thinning.index=thin.index.cs)
		  }
		  else
		  {
		    values <- get.tfr.parameter.traces.cs(mcmc.set$mcmc.list, country.obj, par, 
		                                          burnin=burnin, thinning.index=thin.index)
		  }
			write.values.into.file.cdep(par, values, outdir.thin.mcmc, country.code=country.obj$code,
										compression.type=thinned.mcmc$compression.type)
		}
	}
	if (mcmc.set$meta$nr_countries > mcmc.set$meta$nr_countries_estimation) {
	    cindex <- if(!is.null(update.with.countries)) update.with.countries[update.with.countries > mcmc.set$meta$nr_countries_estimation] else 
	        (mcmc.set$meta$nr_countries_estimation+1):mcmc.set$meta$nr_countries
	    if(length(cindex) > 0)
		    .update.thinned.extras(mcmc.set, cindex, burnin=burnin, nr.points=nr.points, dir=outdir.thin.mcmc, verbose=verbose)
	}
	invisible(structure(list(meta=meta, mcmc.list=list(thinned.mcmc)), class='bayesTFR.mcmc.set'))
}

.update.thinned.extras <- function (mcmc.set, country.index, burnin, nr.points, dir, verbose=TRUE) {
	if(verbose) cat('done.\nStoring country-specific parameters for extra countries ...')
	# thin mcmc for extra countries (they can have different length than the other countries)
	par.names.cs <- tfr.parameter.names.cs(trans=FALSE)
	for (country in country.index){
		if(!(country %in% mcmc.set$meta$id_DL)) next
		country.obj <- get.country.object(country, mcmc.set$meta, index=TRUE)
		for (par in par.names.cs) {
			values <- get.tfr.parameter.traces.cs(mcmc.set$mcmc.list, country.obj, par, 
											burnin=burnin)
			selected.simu <- get.thinning.index(nr.points, dim(values)[1])
			if (length(selected.simu$index) < nr.points)
				selected.simu$index <- sample(selected.simu$index, nr.points, replace=TRUE)
			values <- values[selected.simu$index,]
			write.values.into.file.cdep(par, values, dir, country.code=country.obj$code, 
								compression.type=mcmc.set$meta$compression.type)
		}
	}
	if(verbose) cat('done.\n')
}

.update.meta.for.thinned.mcmc <- function(thin.mcmc.set, new.mcmc.set) {
	# Updating meta object when extra countries are added
	keep.meta <- thin.mcmc.set$meta
	thin.meta <- new.mcmc.set$meta
	for(item in c('output.dir', 'is.thinned', 'parent.iter', 'nr.chains', 'compression.type'))
		thin.meta[[item]] <- keep.meta[[item]]
	thin.meta$parent.meta <- new.mcmc.set$meta
	for (ichain in 1:length(thin.mcmc.set$mcmc.list)) {
		thin.mcmc.set$mcmc.list[[ichain]]$meta <- thin.meta
		store.bayesTFR.object(thin.mcmc.set$mcmc.list[[ichain]], 
				file.path(thin.meta$output.dir, thin.mcmc.set$mcmc.list[[ichain]]$output.dir))
	}
	store.bayesTFR.meta.object(thin.meta, thin.meta$output.dir)
	return(thin.mcmc.set)
}

has.tfr.prediction <- function(mcmc=NULL, sim.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.')
	if(file.exists(file.path(sim.dir, 'predictions', 'prediction.rda'))) return(TRUE)
	return(FALSE)
}

get.tfr.prediction <- function(mcmc=NULL, sim.dir=NULL, mcmc.dir=NULL) {
	############
	# Returns an object of class bayesTFR.prediction
	# Set mcmc.dir to NA, if the prediction object should not have a pointer 
	# to the corresponding mcmc traces.
	############
	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')
	pred.file <- file.path(output.dir, 'prediction.rda')
	if(!file.exists(pred.file)) {
		warning('File ', pred.file, ' does not exist.')
		return(NULL)
	}
	load(file=pred.file)
	bayesTFR.prediction$output.directory <- normalizePath(output.dir)
	
	pred <- bayesTFR.prediction
	if(!is.null(pred$mcmc.set) && is.null(pred$mcmc.set$meta$annual.simulation)) pred$mcmc.set$meta$annual.simulation <- FALSE
	if(!is.null(mcmc.dir) && is.na(mcmc.dir)) return(pred)
	# re-route mcmcs if necessary and load
	if(!is.null(mcmc.dir) || !has.tfr.mcmc(pred$mcmc.set$meta$output.dir)) {
		#if((!is.null(mcmc.dir) && !is.na(mcmc.dir)) || is.null(mcmc.dir)) {
			new.path <- file.path(sim.dir, basename(pred$mcmc.set$meta$output.dir))
			if (has.tfr.mcmc(new.path)) pred$mcmc.set <- get.tfr.mcmc(new.path)
			else {
				est.dir <- if(is.null(mcmc.dir)) sim.dir else mcmc.dir
				pred$mcmc.set <- get.tfr.mcmc(est.dir)
			}
		#}
	}
	return(pred)
}

.get.model.est <- function(excl.what = c(), mcmc.list=NULL, country=NULL, sim.dir=NULL,
                           remove.duplicates = TRUE, country.code = deprecated(), ISO.code = deprecated()) {

    if (lifecycle::is_present(country.code)) {
        lifecycle::deprecate_warn("7.1-0", "tfr.bias.sd(country.code)", "tfr.bias.sd(country)",
                                  details = "The same applies to get.bias.model() and get.sd.model()")
        country <- country.code
    }
    if (lifecycle::is_present(ISO.code)) {
        lifecycle::deprecate_warn("7.1-0", "tfr.bias.sd(ISO.code)", "tfr.bias.sd(country)",
                                  details = "The same applies to get.bias.model() and get.sd.model()")
        country <- ISO.code
    }
    if (is.null(mcmc.list)) 
        mcmc.list <- get.tfr.mcmc(sim.dir)
    if (is.null(mcmc.list)) {
        warning('MCMC does not exist.')
        return(NULL)
    }
    if (is.null(mcmc.list$mcmc.list[[1]]$uncertainty) || !mcmc.list$mcmc.list[[1]]$uncertainty) 
    {
        stop("MCMC does not consider uncertainty of past TFR.")
    }
    if(is.null(country)) stop("Country must be given.")
    meta <- mcmc.list$meta
    country.obj <- get.country.object(country, meta)
    if (country.obj$index %in% meta$extra) model_est <- meta$raw_data_extra[[country.obj$index]]
    else model_est <- meta$raw_data.original[meta$raw_data.original$country_code == country.obj$code,]
    model_est <- model_est[, !(names(model_est) %in% c('country_code', 'Country.or.area', 'year', 'tfr', 'country_index', 'eps', excl.what))]
    if(remove.duplicates)
        model_est <- model_est[!duplicated(model_est),]
    return(list(table = model_est, cntry.index = country.obj$index, meta = meta))
}

get.bias.model <- function(mcmc.list=NULL, country=NULL, sim.dir=NULL, ...) {
    model_est <- .get.model.est(excl.what = "std", mcmc.list = mcmc.list, country = country, sim.dir = sim.dir, ...)
    return(list(model = model_est$meta$bias_model[[model_est$cntry.index]], table=model_est$table))
}

get.std.model <- function(mcmc.list=NULL, country=NULL, sim.dir=NULL, ...) {
    model_est <- .get.model.est(excl.what = "bias", mcmc.list = mcmc.list, country = country, sim.dir = sim.dir, ...)
    return(list(model=model_est$meta$std_model[[model_est$cntry.index]], table=model_est$table))
}

tfr.bias.sd <- function(mcmc.list=NULL, country=NULL, sim.dir=NULL, ...) {
    model_est <- .get.model.est(mcmc.list = mcmc.list, country = country, 
                                sim.dir = sim.dir, remove.duplicates = FALSE, ...)
    covariates <- if (model_est$cntry.index %in% model_est$meta$extra) model_est$meta$extra_covariates[[model_est$cntry.index]] else model_est$meta$covariates

    model_est_dt <- data.table::as.data.table(model_est$table)
    model_est_dt[, 'count':=.N, by=eval(covariates)]
    model_est_df <- as.data.frame(model_est_dt)
    model_est_df <- model_est_df[!duplicated(model_est_df),]
    return(list(model_bias = model_est$meta$bias_model[[model_est$cntry.index]], 
                model_sd = model_est$meta$std_model[[model_est$cntry.index]], 
                table = model_est_df))
}

has.est.uncertainty <- function(mcmc.set)
    return(length(mcmc.set$mcmc.list)>0 && !is.null(mcmc.set$mcmc.list[[1]]$uncertainty) && mcmc.set$mcmc.list[[1]]$uncertainty)

get.tfr.estimation <- function(mcmc.list=NULL, country=NULL, sim.dir=NULL, 
                               burnin=0, thin = 1, probs=NULL, adjust=TRUE, 
                               country.code = deprecated(), ISO.code = deprecated()) {

    if (lifecycle::is_present(country.code)) {
        lifecycle::deprecate_warn("7.1-0", "get.tfr.estimation(country.code)", "get.tfr.estimation(country)")
        country <- country.code
    }
    if (lifecycle::is_present(ISO.code)) {
        lifecycle::deprecate_warn("7.1-0", "get.tfr.estimation(ISO.code)", "get.tfr.estimation(country)")
        country <- ISO.code
    }
    
    if(is.character(mcmc.list)) {
        sim.dir <- mcmc.list
        mcmc.list <- NULL
    }
  if (is.null(mcmc.list)) 
    mcmc.list <- get.tfr.mcmc(sim.dir)
  if (is.null(mcmc.list)) {
    warning('MCMC does not exist.')
    return(NULL)
  }
    if (!has.est.uncertainty(mcmc.list))
    {
        stop("MCMC does not consider uncertainty of past TFR.")
    }
    if (is.null(country))
        stop("Country  must be given.")

  country.obj <- get.country.object(country, mcmc.list$meta)
  tfr_table <- get.tfr.parameter.traces.cs(mcmc.list$mcmc.list, country.obj, 'tfr', burnin = burnin, thin=thin)
  if (adjust)
  {
    shift <- get.tfr.shift.estimation(country.obj$code, mcmc.list$meta)
    if (!is.null(shift)) tfr_table <- t(t(tfr_table) + shift)
  }
  output <- list(tfr_table=tfr_table, country.obj = country.obj)
  if (!is.null(probs))
  {
    tfr_quantile <- apply(tfr_table, 2, quantile, probs = probs)
    if (is.null(dim(tfr_quantile))) tfr_quantile <- matrix(tfr_quantile, nrow = 1)
    tfr_quantile <- data.table::data.table(t(tfr_quantile))
    start.year <- mcmc.list$meta$start.year
    end.year <- mcmc.list$meta$present.year
    if (mcmc.list$meta$annual.simulation)
      tfr_quantile$year <- start.year:end.year
    else
      tfr_quantile$year <- seq(start.year+3, end.year-2, 5)
    output$tfr_quantile <- tfr_quantile
  }
  class(output) <- "bayesTFR.estimation"
  return(output)
}

print.bayesTFR.estimation <- function(x, ...){
    res <- list(country.obj = x$country.obj, tfr_table = head(x$tfr_table))
    if(!is.null(x$tfr_quantile)) res$tfr_quantile <- head(x$tfr_quantile)
    print(res)
}

get.regtfr.prediction <- function(sim.dir, country=NULL){
  dir <- file.path(sim.dir, 'subnat')
  if(!file.exists(dir)) stop("No subnational predictions in ", sim.dir)
  cdirs <- list.files(dir, pattern='^c[0-9]+', full.names=FALSE)
  if(length(cdirs)==0) stop("No subnational predictions in ", sim.dir)
  if (!is.null(country)) {
    cdirs <- cdirs[cdirs == paste0("c", country)]
    if(length(cdirs)==0) stop("No subnational predictions for country ", country, " in ", sim.dir)
    return(get.tfr.prediction(sim.dir=file.path(dir, cdirs[1]), mcmc.dir=NA))
  }
  preds <- NULL
  for (d in cdirs) 
    preds[[sub("c", "", d)]] <- get.tfr.prediction(sim.dir=file.path(dir, d), mcmc.dir=NA)
  return(preds)
}

.do.get.convergence.all <- function(type, package, sim.dir) {
	diag.dir <- file.path(sim.dir, 'diagnostics')
	if (!file.exists(diag.dir)) return(NULL)
	files <- list.files(diag.dir, pattern=paste('^', package, '[.]convergence_[0-9]+_[0-9]+[.]rda$', sep=''), 
							full.names=FALSE)
	if(length(files) <= 0) return(NULL)
	thin.matches <- regexpr('e_[0-9]+', files)
	bi.matches <- regexpr('[0-9]+[.]', files)
	result <- list()
	for(i in 1:length(files)) {
		thin <- as.integer(substr(files[i], start=thin.matches[i]+2, 
								stop=thin.matches[i]+attr(thin.matches, 'match.length')[i]-1))
		burnin <- as.integer(substr(files[i], start=bi.matches[i], 
								stop=bi.matches[i]+attr(bi.matches, 'match.length')[i]-2))
		result[[i]] <- do.call(paste('get.', type, '.convergence', sep=''), 
								list(sim.dir, thin=thin, burnin=burnin))
	}
	return(result)
}
get.tfr.convergence.all <- function(sim.dir=file.path(getwd(), 'bayesTFR.output')) {
	return(.do.get.convergence.all('tfr', 'bayesTFR', sim.dir=sim.dir))
}

get.tfr.convergence <- function(sim.dir=file.path(getwd(), 'bayesTFR.output'), 
									thin=80, burnin=2000) {
	file.name <- file.path(sim.dir, 'diagnostics', paste('bayesTFR.convergence_', thin, '_', burnin, '.rda', sep=''))
	if(!file.exists(file.name)){
		warning('Convergence diagnostics in ', sim.dir, ' for burnin=', burnin, 
					' and thin=', thin, ' does not exist.')
		return(NULL)
	}
	bayesTFR.convergence <- local({
								load(file.name)
								bayesTFR.convergence})
	return(bayesTFR.convergence)
}

get.tfr3.convergence.all <- function(sim.dir=file.path(getwd(), 'bayesTFR.output')) {
	return(get.tfr.convergence.all(sim.dir=file.path(sim.dir, 'phaseIII')))
}

get.tfr3.convergence <- function(sim.dir=file.path(getwd(), 'bayesTFR.output'), 
									thin=60, burnin=10000) {
	return(get.tfr.convergence(file.path(sim.dir, 'phaseIII'), thin=thin, burnin=burnin))
}



has.mcmc.converged <- function(diag) return(diag$status['green'])	
get.burned.tfr.traces <- function(mcmc, par.names, burnin=0, thinning.index=NULL) {
	# get traces that are already loaded in the mcmc object
	traces <- mcmc$traces[, par.names, drop = FALSE]
	discard <- burnin - mcmc$traces.burnin
	if (discard > 0)
		traces <- traces[-seq(1, discard),,drop=FALSE]
	if(!is.null(thinning.index))
		traces <- traces[thinning.index,,drop=FALSE]
	return(traces)
}
"bdem.parameter.traces" <- function(mcmc, ...) UseMethod("bdem.parameter.traces")

bdem.parameter.traces.bayesTFR.mcmc <- function(mcmc, par.names, ...) {
	tran.names <- totran.names <- all.standard.names <- backtran.names <- tobacktran.names <- c()
	# Load traces from the disk
	if(is.null(mcmc$meta$phase) || mcmc$meta$phase == 2) {
		all.standard.names <- c(tfr.parameter.names(meta=mcmc$meta), get.trans.parameter.names(), 
							tfr.parameter.names.cs(), get.trans.parameter.names(cs=TRUE), 'tfr')
		tran.names <- c(get.trans.parameter.names(), get.trans.parameter.names(cs=TRUE))
		totran.names <- c(get.totrans.parameter.names(), get.totrans.parameter.names(cs=TRUE))
		backtran.names <- get.backtrans.parameter.names(cs=TRUE)
		tobacktran.names <- get.tobacktrans.parameter.names(cs=TRUE)
	} else {
		if(mcmc$meta$phase == 3) all.standard.names <- c(tfr3.parameter.names(), tfr3.parameter.names.cs())
	}
	return(.do.get.traces(mcmc, par.names=par.names, ..., 
							all.standard.names=all.standard.names,
							tran.names=tran.names, totran.names=totran.names,
							backtran.names=backtran.names, tobacktran.names=tobacktran.names))
}

.do.get.traces <- function(mcmc, par.names, file.postfix='', par.names.postfix='', burnin=0, 
							thinning.index=NULL, all.standard.names=c(), tran.names=c(), totran.names=c(), 
							backtran.names=c(), tobacktran.names=c()) {
  if (length(par.names) == 0) return (NULL)
  tran.names.l <- nchar(tran.names)
	ltran.names <- length(tran.names)
	totran.names.l <- nchar(totran.names)
	has.tran <- rep(FALSE, ltran.names)
	lbacktran.names <- length(backtran.names)
	backtran.names.l <- nchar(backtran.names)
	has.backtran <- rep(FALSE, lbacktran.names)
	values <- c()
 	valnames <- c()
 	loaded.files <- c()
 	if(is.null(mcmc$meta$package.version) || mcmc$meta$package.version < "7.0.0") 
 	    compr.settings <- .get.compression.settings.obsolete(mcmc$compression.type)
 	else
 	    compr.settings <- .get.compression.settings(mcmc$compression.type)
 	for(name in par.names) {
 		if(lbacktran.names > 0) {
 			new.has.backtran <- backtran.names %in% sapply(backtran.names.l, function(x) substr(name, 1, x))
 			has.backtran <- has.backtran | new.has.backtran
 			if(any(new.has.backtran)) next  # name is a back-transformation variable, therefore skip to next one
 		} 
 		if(ltran.names > 0) {
 			new.has.tran <- tran.names %in% sapply(tran.names.l, function(x) substr(name, 1, x))
 			has.tran <- has.tran | new.has.tran
 			if(any(new.has.tran)) next  # name is a transformation variable, therefore skip to next one
 		}
 		if (any(name == valnames)) next # name already loaded
 			
 		if (!any(name==all.standard.names)) {
 			#trim the name to the standard names, e.g. make 'gamma' from 'gamma_3'
 			name.to.load <- NULL
 			for (sname in all.standard.names) {
 				if (length(grep(paste('^', sname, '_', sep=''), name)) > 0) {
 					name.to.load <- sname
 					break
 				}
 			}
 			if (is.null(name.to.load)) {
 				warning('Parameter ', name, ' not found.')
 				next
 			}
 		} else {name.to.load <- name}
 		if (any(name.to.load == valnames)) next # name already loaded
 		#file.name <- paste(name.to.load, file.postfix, '.txt.bz2',sep='')
 		file.name <- paste(name.to.load, file.postfix, '.txt', compr.settings[2], sep='') 		
 		if (any(file.name == loaded.files)) next
 		if(is.null(mcmc$meta$package.version) || mcmc$meta$package.version < "7.0.0" || 
 		   compr.settings[1] %in% c("xzfile", "bzfile")) { # old ways of reading files
 		    #con <- bzfile(file.path(mcmc$meta$output.dir, mcmc$output.dir, file.name), open="rb")
 		    #con <- file(file.path(mcmc$meta$output.dir, mcmc$output.dir, file.name), open="r")
 		    con <- do.call(compr.settings[1], list(file.path(mcmc$meta$output.dir, mcmc$output.dir, file.name), 
 								open=paste("r", compr.settings[3], sep='')))
 		    if(compr.settings[3]=='b')  # binary connection
 		    	raw.con <- textConnection(readLines(con))
 		    else raw.con <- con
 		    # close(con)
		    vals <- as.matrix(read.table(file = raw.con))
		    close(raw.con)
		    #vals <- as.matrix(read.table(file = con))
		    if(compr.settings[3]=='b') close(con)
 		} else # new way of reading files using data.table
 		    vals <- as.matrix(data.table::fread(file.path(mcmc$meta$output.dir, mcmc$output.dir, file.name)))
		loaded.files <- c(loaded.files, file.name)
		if (dim(vals)[2] > 1) { #2d parameters get postfix
			valnames <- c(valnames, paste(name.to.load, '_', 1:dim(vals)[2], par.names.postfix, sep=''))
		} else {
			valnames <- c(valnames, paste(name.to.load, par.names.postfix, sep=''))
		}
		if (burnin > 0) {
			if (burnin > dim(vals)[1]) stop('Burnin is larger than the data size.')
			vals <- vals[-seq(1, burnin),,drop=FALSE]
		}
		if(!is.null(thinning.index))
			vals <- vals[thinning.index,,drop=FALSE]
		values <- cbind(values, vals)
	} # end of loading loop
	#stop('')
 	colnames(values) <- valnames
	if(ltran.names == 0 && lbacktran.names == 0 ) return(values)
	# get transformed variables (alpha, delta, gamma)
	if(ltran.names > 0 && any(has.tran)) {
		for (i in which(has.tran)) {
			full.names <- grep(paste0(totran.names[i],'_'), valnames, value=TRUE)
			if (length(grep(paste0('^', totran.names[i], '(_|$)'), valnames)) == 0) { # load alpha/delta/gamma
				vals <- bdem.parameter.traces(mcmc, c(totran.names[i]), 
								file.postfix=file.postfix, par.names.postfix=par.names.postfix,
								burnin=burnin, thinning.index=thinning.index)
				full.names <- grep(paste0(totran.names[i],'_'), c(valnames, colnames(vals)), value=TRUE)
			} else {
				vals <- values[, full.names, drop=FALSE]
			}
			vals <- exp(vals)/apply(exp(vals),1,sum)
			values <- cbind(values, vals)
			valnames <- c(valnames, paste0(tran.names[i],'_', 
											substr(full.names, totran.names.l[i]+2, nchar(full.names))))
			colnames(values) <- valnames
		}
	}
	#stop('')
	# get back-transformed variables (Triangle_c1-3)
	if(lbacktran.names > 0 && any(has.backtran)) {
	  vals <- c()
		for(i in 1:length(tobacktran.names)) {		
			full.names <- grep(paste0(tobacktran.names[i],'_'), valnames, value=TRUE)
			if (length(grep(paste0('^', tobacktran.names[i], '(_|$)'), valnames)) == 0) { # load gamma, U, Triangle_c4
				vals <- cbind(vals, bdem.parameter.traces(mcmc, c(tobacktran.names[i]), 
								file.postfix=file.postfix, par.names.postfix=par.names.postfix,
								burnin=burnin, thinning.index=thinning.index))
				full.names <- grep(paste0(tobacktran.names[i],'_'), c(valnames, colnames(vals)), value=TRUE)
			} else {
				vals <- cbind(vals, values[, full.names, drop=FALSE])
			}
		}
		vals <- .get.backtransformed.Triangles(vals, backtran.names[has.backtran])
		values <- cbind(values, vals)
	}
	return(values)
}

.get.backtransformed.Triangles <- function(data, par.names) {
	data.par.names <- colnames(data)
	# set names of the parameter columns for this country
	U.var <- grep('^U_c', data.par.names, value=TRUE) 
	d.var <- grep('^d_c', data.par.names, value=TRUE) 
	Triangle_c4.var <- grep("^Triangle_c4_c", data.par.names, value=TRUE) 
	gamma.vars <- sapply(paste0('^gammat_',1:3), function(x) grep(x, data.par.names, value=TRUE))
	# transform gamma_ci into Triangle_ci 
	# compute p_ci 
	#p_ci <- matrix(NA, nrow(data), 3) 
	#p_ci_denominator <- apply(exp(data[,gamma.vars]), 1, sum) 
	#for (i in 1:3) p_ci[,i] <- exp(data[,gamma.vars[i]])/p_ci_denominator 
	Delta <- data.frame(
    	Triangle_c1 = (data[,U.var] - data[,Triangle_c4.var])*data[,gamma.vars[1]],
    	Triangle_c2 = (data[,U.var] - data[,Triangle_c4.var])*data[,gamma.vars[2]],
    	Triangle_c3 = (data[,U.var] - data[,Triangle_c4.var])*data[,gamma.vars[3]]
    ) 
    Delta <- Delta[,par.names, drop=FALSE]
    suffix <- strsplit(data.par.names, "_")[[1]]
    suffix <- suffix[length(suffix)]
    colnames(Delta) <- paste(colnames(Delta), suffix, sep="_")
	return(as.matrix(Delta))
	
}


get.all.parameter.names <- function() {
    # First element in each tuple indicates if the parameter is transformable, 
    # the second how many values it has.
	return(list(alpha=c(TRUE, 3),  delta=c(FALSE, 3), Triangle4=c(FALSE, 1), delta4=c(FALSE, 1), 
				psi=c(FALSE, 1), chi=c(FALSE, 1), a_sd=c(FALSE, 1), b_sd=c(FALSE, 1), 
				const_sd=c(FALSE, 1), S_sd=c(FALSE, 1), sigma0=c(FALSE, 1), mean_eps_tau=c(FALSE, 1),
				sd_eps_tau=c(FALSE, 1)))
}				

get.all.parameter.names.cs <- function() {
    # First element in each tuple indicates if the parameter is transformable, 
    # the second how many values it has.
	return(list(gamma=c(1, 3), U=c(0, 1), d=c(0, 1), Triangle_c4=c(0, 1), 
				Triangle_c1=c(2, 1), Triangle_c2=c(2, 1), Triangle_c3=c(2, 1)))
}

get.totrans.parameter.names <- function(cs=FALSE) {
	pars <- if(cs) get.all.parameter.names.cs() else get.all.parameter.names() 
	is.trans <- sapply(pars, function(x) return(x[1] == 1))
	return(names(pars)[is.trans])
}

get.trans.parameter.names <- function(cs=FALSE) return(paste(get.totrans.parameter.names(cs), 't', sep=''))

get.backtrans.parameter.names <- function(cs=FALSE) {
	if(!cs) return (c())
	pars <- get.all.parameter.names.cs()
	is.backtrans <- sapply(pars, function(x) return(x[1] == 2))
	return(names(pars)[is.backtrans])
}

get.tobacktrans.parameter.names <- function(cs=FALSE) {
	if(!cs) return (c())
	return(c("U", "gammat", "Triangle_c4"))
}

get.other.parameter.names <- function(cs=FALSE) {
	pars <- if(cs) get.all.parameter.names.cs() else get.all.parameter.names()
	is.trans <- sapply(pars, function(x) return(x[1] > 0))
	return(names(pars)[!is.trans])
}

get.all.parameter.names.extended <- function(cs=FALSE) {
  pars <- c()
	all.pars <- if(cs) get.all.parameter.names.cs() else get.all.parameter.names()
	for (ipar in 1:length(all.pars)) {
		par <- all.pars[ipar]
		paropt <- unlist(par)
		name <- names(par)
		pars <- c(pars, if (paropt[2] > 1) paste(name, 1:paropt[2], sep='_') else name)
		if (paropt[1]) {
			name <- paste(name,'t', sep='')
			pars <- c(pars, if (paropt[2] > 1) paste(name, 1:paropt[2], sep='_') else name)
		}
	}
	if (cs) pars <- c(pars, 'tfr')
	else pars <- c(pars, "rho_phase2")
	return(pars)
}

tfr.parameter.names.extended <- function() {
	# return a list of all parameters with those parameters extended by i
	# that have more than 1 value, e.g. alpha_2, alphat_1, delta_3
	return(get.all.parameter.names.extended(cs=FALSE))
}

tfr.parameter.names.cs.extended <- function(country.code=NULL) {
	# return a list of cs-parameters with gamma extended by i
	pars <- get.all.parameter.names.extended(cs=TRUE)
	if(!is.null(country.code)) pars <- paste(pars, '_c', country.code, sep='')
	return(pars)
}

tfr.parameter.names <- function(trans=NULL, meta=NULL) {
	# Return all country-independent parameter names. 
	# trans can be NULL or logical.
	# If 'trans' is TRUE,
	# names of the transformable parameters (alpha) are replaced by the names 
	# of the transformed parameters.
	# If 'trans' is NULL, all parameter names, 
	# including the transformable parameters are returned.
	other.pars <- get.other.parameter.names(cs=FALSE)
	if (!is.null(meta) && !is.null(meta$ar.phase2) && meta$ar.phase2) other.pars <- c(other.pars, 'rho_phase2')
	if (is.null(trans)) return (c(get.totrans.parameter.names(), get.trans.parameter.names(), other.pars))
	if (trans) return (c(get.trans.parameter.names(), other.pars))
	return (c(get.totrans.parameter.names(), other.pars))
}
	
tfr.parameter.names.cs <- function(country.code=NULL, trans=NULL, back.trans=TRUE) {
	#Return all country-specific parameter names. 
	# See comments in tfr.parameter.names(). Transformable parameter is gamma.
	#If country is not NULL, it must be a country code.
	#It is attached to the parameter name.
	par.names <- get.other.parameter.names(cs=TRUE)
	backtrans.names <- if(back.trans) get.backtrans.parameter.names(cs=TRUE) else c()
	if (is.null(trans)) par.names <- c(par.names, backtrans.names,
										get.totrans.parameter.names(cs=TRUE), 										
										get.trans.parameter.names(cs=TRUE))
	else par.names <- c(par.names, backtrans.names,
								if (trans) get.trans.parameter.names(cs=TRUE) 
									else get.totrans.parameter.names(cs=TRUE))
	if (is.null(country.code))
		return(par.names)
	return(paste(par.names, '_c', country.code, sep=''))
}

tfr3.parameter.names <- function() return(c('mu', 'rho', 'sigma.mu', 'sigma.rho', 'sigma.eps'))
tfr3.parameter.names.cs <- function(country.code=NULL) return(paste(c('mu', 'rho'), '.c', country.code,sep=''))

get.total.iterations <- function(mcmc.list, burnin=0) {
	# Return total number of iterations sum-up over chains after discarding burnin in each chain
	get.iter <- function(x) return(x$finished.iter - burnin)
	return(sum(sapply(mcmc.list, get.iter)))
}

get.thinned.burnin <- function(mcmc, burnin) {
	if (burnin==0) return(0)
	if (mcmc$thin == 1) return(burnin)
	return(1 + if(burnin >= mcmc$thin) length(seq(mcmc$thin, burnin, by=mcmc$thin)) else 0)
}

get.stored.mcmc.length <- function(mcmc.list, burnin=0) {
	# Return total number of iterations sum-up over chains after discarding burnin in each chain,
	# taking into account the original value of thin.
	# It should correspond to the total length of all chains stored on disk, minus burnin.
	get.iter <- function(x) return(x$length - get.thinned.burnin(x, burnin))
	return(sum(sapply(mcmc.list, get.iter)))
}

get.stored.mcmc.length.extra <- function(meta, countries.index = NULL, nr.chains = 1, burnin=0) {
    if (is.null(meta$extra) || is.null(countries.index) || !any(countries.index %in% meta$extra)) return(NA)
    stored.iter <- NA
    for (country in countries.index){
        if(! country %in% meta$extra) next
        stored.iter <- min(stored.iter, 
                           floor((meta$extra_iter[country] - burnin) / meta$extra_thin[country]) * nr.chains, na.rm = TRUE)
    }
    return(stored.iter)
}


do.get.tfr.parameter.traces <- function(is.cs, mcmc.list, par.names, country.obj=NULL, 
										burnin=0, thinning.index=NULL, thin=NULL) {
	# get parameter traces either from disk or from memory (if they were already loaded)
	# par.names are either country-independent (if is.cs is FALSE), or country-specific (if is.cs is TRUE)
  values <- c()
	if (is.null(thinning.index) && !is.null(thin) && thin > mcmc.list[[1]]$thin) {
		total.iter <- get.stored.mcmc.length(mcmc.list, burnin)
		thinning.index <- unique(round(seq(1, total.iter, by=thin/mcmc.list[[1]]$thin)))
	}
	at.iter <- 1
	for(mcmc in mcmc.list) {
		this.thinning.index <- NULL
		if (!is.null(country.obj) && !is.null(mcmc$meta$extra) && country.obj$index %in% mcmc$meta$extra)
		{
		  th.burnin <- ceiling(burnin / mcmc$meta$extra_thin[country.obj$index])
		}
		else
		{
		  th.burnin <- get.thinned.burnin(mcmc, burnin)
		}
		if(!is.null(thinning.index)) {
		  step.size <- mcmc$length
		  if (!is.null(country.obj) && !is.null(mcmc$meta$extra) && country.obj$index %in% mcmc$meta$extra)
		  {
		    step.size <- floor(mcmc$meta$extra_iter[country.obj$index] /mcmc$meta$extra_thin[country.obj$index])
		  }
			this.thinning.index <- thinning.index[(thinning.index >= at.iter) & 
									(thinning.index < at.iter+step.size-th.burnin)] - at.iter+1
			if (length(this.thinning.index) == 0) {
				at.iter <- at.iter+step.size-th.burnin
				next
			}
		}
    	if (no.traces.loaded(mcmc) || th.burnin < mcmc$traces.burnin) {
    		traces <- if(is.cs) load.tfr.parameter.traces.cs(mcmc, country.obj$code, par.names, burnin=th.burnin, 
    											thinning.index=this.thinning.index)
    					else bdem.parameter.traces(mcmc, par.names, burnin=th.burnin, thinning.index=this.thinning.index)
        } else {
            traces <- if(is.cs) get.burned.tfr.traces(mcmc, get.full.par.names.cs(par.names, colnames(mcmc$traces), 
            									country=country.obj$index), 
            								th.burnin, thinning.index=this.thinning.index)
            		  else get.burned.tfr.traces(mcmc, par.names, th.burnin, thinning.index=this.thinning.index)
        }
       	values <- rbind(values, traces)
       	if(!is.null(thinning.index))
       	{
       	  at.iter <- at.iter + step.size - th.burnin
       	}
       	else at.iter <- at.iter+mcmc$length-th.burnin
    }
    return(values)
}

get.tfr.parameter.traces <- function(mcmc.list, par.names=tfr.parameter.names(), 
										burnin=0, thinning.index=NULL, thin=NULL) {
	# get parameter traces either from disk or from memory, if they were already loaded
	return(do.get.tfr.parameter.traces(is.cs=FALSE, mcmc.list=mcmc.list, par.names=par.names, 
										burnin=burnin, thinning.index=thinning.index, thin=thin))
}

get.tfr.parameter.traces.cs <- function(mcmc.list, country.obj, par.names=tfr.parameter.names.cs(), 
										burnin=0, thinning.index=NULL, thin=NULL) {
	# country.obj is result of get.country.object()
	# get traces for country-specific parameters either from disk or from memory, if they were already loaded
	return(do.get.tfr.parameter.traces(is.cs=TRUE, mcmc.list=mcmc.list, par.names=par.names, country.obj=country.obj,
										burnin=burnin, thinning.index=thinning.index, thin=thin))
}

get.tfr3.parameter.traces <- function(mcmc.list, par.names=tfr3.parameter.names(), ...)
  return(get.tfr.parameter.traces(mcmc.list, par.names, ...))
	
get.tfr3.parameter.traces.cs <- function(mcmc.list, country.obj, par.names=tfr3.parameter.names.cs(), ...)
	return(get.tfr.parameter.traces.cs(mcmc.list, country.obj, par.names, ...))

load.tfr.parameter.traces <- function(mcmc, par.names=tfr.parameter.names(), burnin=0, thinning.index=NULL) 
 	return(bdem.parameter.traces(mcmc, par.names, burnin=burnin, thinning.index=thinning.index))

load.tfr.parameter.traces.cs <- function(mcmc, country, par.names=tfr.parameter.names.cs(), burnin=0, 
										thinning.index=NULL) {
  return(bdem.parameter.traces(mcmc, par.names, paste('_country', country, sep=''),
						par.names.postfix=paste('_c', country, sep=''), burnin=burnin, 
						thinning.index=thinning.index))
}

load.tfr3.parameter.traces <- function(mcmc, par.names=tfr3.parameter.names(), ...)
	return(load.tfr.parameter.traces(mcmc, par.names=par.names, ...))

load.tfr3.parameter.traces.cs <- function(mcmc, country, par.names=tfr3.parameter.names.cs(), ...)
	return(load.tfr.parameter.traces.cs(mcmc, country=country, par.names=par.names, ...))

load.tfr.parameter.traces.all <- function(mcmc, par.names=tfr.parameter.names(), 
										 par.names.cs=tfr.parameter.names.cs(),
										 burnin=0, thinning.index=NULL) {
	result <- load.tfr.parameter.traces(mcmc, par.names, burnin=burnin, thinning.index=thinning.index)
	if(!is.null(par.names.cs))
		for (country in get.countries.index(mcmc$meta)) {
			result <- cbind(result, 
						load.tfr.parameter.traces.cs(mcmc, 
												    get.country.object(country, 
												         mcmc$meta, index=TRUE)$code, 
												    par.names.cs, burnin=burnin,
												    thinning.index=thinning.index))
		}
	return (result)
}

load.tfr3.parameter.traces.all <- function(mcmc, par.names=tfr3.parameter.names(), 
										 par.names.cs=tfr3.parameter.names.cs(), ...)
	return(load.tfr.parameter.traces.all(mcmc, par.names=par.names, par.names.cs=par.names.cs, ...))
	
get.full.par.names.cs <- function(par.names, full.par.names, country=NULL, index=FALSE) {
	# Return full name of par.names that are included in full.par.names
	# which are suppose to be all country-specific parameters.
	# E.g. for 'gamma', it would return 'gamma_1_c1', ..., gamma_1_cX'
	# If index is TRUE, return the index of the matches.
	result <- c()	
	for (name in par.names) {
		pattern <- paste('^', name, '_.*c', sep='')
		if (!is.null(country)) {
			pattern <- paste(pattern, country, '$', sep='')
		}
		result <- c(result, grep(pattern, full.par.names, value=!index))
	}
	return(result)
}

get.full.par.names <- function(par.names, full.par.names, index=FALSE) {
	# Return full name of par.names that are included in full.par.names
	# which are suppose to be all country-independent parameters.
	# E.g. for 'alpha', it would return 'alpha_1', alpha_2, alpha_3'
	# If index is TRUE, return the index of the matches.
	result <- c()	
	for (name in par.names) {
		pattern <- paste('^', name, '(_|$)', sep='') # either has '_' suffix or matches exactly
		result <- c(result, grep(pattern, full.par.names, value=!index))
	}
	return(result)
}

"coda.mcmc" <- function(mcmc, ...) UseMethod("coda.mcmc")

coda.mcmc.bayesTFR.mcmc <- function(mcmc, country=NULL, par.names=NULL, 
						par.names.cs=NULL, burnin=0, thin=1, ...
						) {
	# Return a coda object for this mcmc and parameter names
	index <- NULL
	btobject <- burn.and.thin(mcmc, burnin=burnin, thin=thin)
	thin <- btobject$thin
	th.burnin <- btobject$burnin
	if(!is.null(btobject$index)) index <- btobject$index
	if(is.null(country) && missing(par.names)) 
		par.names <- if(!is.null(mcmc$meta$phase) && mcmc$meta$phase == 3)
			tfr3.parameter.names() else tfr.parameter.names(meta = mcmc$meta)
	if(missing(par.names.cs)) 
		par.names.cs <- if(!is.null(mcmc$meta$phase) && mcmc$meta$phase == 3)
			tfr3.parameter.names.cs() else tfr.parameter.names.cs()

	if (!is.null(country)) { # for specific country
		if (burnin < mcmc$traces.burnin || no.traces.loaded(mcmc)) {
			values <- load.tfr.parameter.traces(mcmc, par.names, burnin=th.burnin, thinning.index=index)
			values <- cbind(values, load.tfr.parameter.traces.cs(mcmc, 
								get.country.object(country, mcmc$meta)$code, par.names.cs, 
								burnin=th.burnin, thinning.index=index))
		} else {
			values <- get.burned.tfr.traces(mcmc, 
							c(get.full.par.names(par.names, colnames(mcmc$traces)), 
							  get.full.par.names.cs(par.names.cs, colnames(mcmc$traces), 
							  	get.country.object(country, mcmc$meta)$code)), th.burnin,
							  	thinning.index=index)
		}
	} else { #no country specified
		if (no.traces.loaded(mcmc)) { # traces not loaded
			values <- load.tfr.parameter.traces.all(mcmc, par.names, par.names.cs, burnin=th.burnin, 
												    thinning.index=index)
		} else { # traces loaded but get the right burnin
			values <- get.burned.tfr.traces(mcmc, c(get.full.par.names(par.names, colnames(mcmc$traces)), 
                             get.full.par.names.cs(par.names.cs, colnames(mcmc$traces))), th.burnin,
                             thinning.index=index)
        }
    }
    # filter out unnecessary parameters
    values <- filter.traces(values, c(par.names, par.names.cs))
    return(mcmc(values, end=mcmc$finished.iter, thin=thin, ...))
}
	
	
coda.list.mcmc <- function(mcmc=NULL, country=NULL, chain.ids=NULL,
							sim.dir=file.path(getwd(), 'bayesTFR.output'), 
							par.names=tfr.parameter.names(), 
							par.names.cs=tfr.parameter.names.cs(), 
							rm.const.pars=FALSE,
							burnin=0, 
							low.memory=FALSE, ...) {
	# return a list of mcmc objects that can be analyzed using the coda package
	mcmc.list <- mcmc
	if (is.null(mcmc.list)) {
		mcmc.list <- stop.if.mcmc.missing(get.tfr.mcmc(sim.dir, chain.ids=chain.ids, low.memory=low.memory, 
									burnin=burnin)$mcmc.list)
	} else {
		mcmc.list <- get.mcmc.list(mcmc.list)
		if (!is.null(chain.ids)) {
			mcmc.list <- mcmc.list[chain.ids]
		}
	}
	result <- list()
	i <- 1
	for(mc in mcmc.list) {
		result[[i]] <- coda.mcmc(mc, country=country, par.names=par.names, par.names.cs=par.names.cs,
								 burnin=burnin, ...)

		if (rm.const.pars) {
			# remove parameters that are constant for all iterations
			if (is.null(dim(result[[i]]))) { # one parameter in the chain
				if(all(result[[i]]==result[[i]][1]))  { # no parameters are kept
					result[[i]] <- NULL
					i <- i - 1
				}
			} else { # multiple parameters in the chain
				if (dim(result[[i]])[2]==1) {
					colindex <- !all(result[[i]][,1] == result[[i]][1,1])
				} else {
					colindex <- !apply(t(apply(result[[i]], 1, '==', result[[i]][1,])), 2, all)
				}
				if (sum(colindex) == 0) { # no parameters are kept
					result[[i]] <- NULL
					i <- i - 1
				} else {
					result[[i]] <- result[[i]][,colindex]
				}
			}
		}
		i <- i+1
	}
	return(mcmc.list(result))
}

coda.list.mcmc3 <- function(mcmc=NULL, country=NULL, chain.ids=NULL,
							sim.dir=file.path(getwd(), 'bayesTFR.output'), 
							par.names=tfr3.parameter.names(), 
							par.names.cs=tfr3.parameter.names.cs(), 
							burnin=0, low.memory=FALSE, ...) {
	if (is.null(mcmc)) 
		mcmc <- get.tfr3.mcmc(sim.dir, chain.ids=chain.ids, low.memory=low.memory, 
									burnin=burnin)$mcmc.list
	else
		if(inherits(mcmc, 'bayesTFR.prediction'))
			stop('Function not available for bayesTFR.prediction objects.')
	return(coda.list.mcmc(mcmc=mcmc, country=country, chain.ids=chain.ids, sim.dir=NULL, 
			par.names=par.names, par.names.cs=par.names.cs, rm.const.pars=FALSE, burnin=burnin, ...))										
}

filter.traces <- function(values, par.names) {
	valuenames <- colnames(values)
    lpar.names <- length(par.names)
    include <- rep(FALSE, length(valuenames))
    for(i in 1:lpar.names) {
    	idx <- grep(paste('^', par.names[i], '$', sep=''), valuenames) # exact match
    	if (length(idx) > 0) {
    		include[idx] <- TRUE
    		next
    	}
    	idx <- grep(paste0('^', par.names[i], '_'), valuenames) # partial match
    	if (length(idx) > 0) include[idx] <- TRUE
    }
	#v <- array(values[,include], c(nrow(values), sum(include)))
	v <- values[,which(include), drop=FALSE]
	#colnames(v) <- colnames(values)[include]
    return(v)
}

stop.if.mcmc.missing <- function(x, ignore.missing = FALSE) {
  if(!ignore.missing && length(x) == 0) stop("No MCMCs found")
  return(x)
}

"get.mcmc.list" <- function(mcmc.list, ...) UseMethod("get.mcmc.list")

get.mcmc.list.bayesTFR.mcmc.set <- function(mcmc.list, ...) return(stop.if.mcmc.missing(mcmc.list$mcmc.list, ...))
get.mcmc.list.bayesTFR.mcmc <- function(mcmc.list, ...) return(stop.if.mcmc.missing(list(mcmc.list), ...))
get.mcmc.list.bayesTFR.prediction <- function(mcmc.list, ...) return(stop.if.mcmc.missing(mcmc.list$mcmc.set$mcmc.list, ...))
get.mcmc.list.list <- function(mcmc.list, ...) return(stop.if.mcmc.missing(mcmc.list, ...))

"get.mcmc.meta" <- function(meta, ...) UseMethod("get.mcmc.meta")
get.mcmc.meta.bayesTFR.mcmc.set <- function(meta, ...) return(meta$meta)
get.mcmc.meta.bayesTFR.mcmc.meta <- function(meta, ...) return(meta)
get.mcmc.meta.bayesTFR.mcmc <- function(meta, ...) return(meta$meta)
get.mcmc.meta.bayesTFR.prediction <- function(meta, ...) return(meta$mcmc.set$meta)


get.country.object <- function(country, meta=NULL, country.table=NULL, index=FALSE) {
	# If meta object is not given, country.table containing columns 'code' and 'name' must be given. 
	# If 'country' is numeric, 'index' determines if 'country' is an index (TRUE) or code (FALSE)
	if (!is.null(meta)) {
		codes <- meta$regions$country_code
		names <- meta$regions$country_name
	} else {
		codes <- country.table[,'code']
		names <- country.table[,'name']	
	}
	l <- length(codes)
	found <- TRUE
	if (is.numeric(country)) {
		if (index) { # index
			country.idx <- country
			country.code <- codes[country.idx]
		} else { # country code
			country.code <- country
			country.idx <- (1:l)[codes==country.code]
		}
		country.name <- as.character(names[country.idx])
		if (length(country.name) == 0) found <- FALSE
	} else { # country is character string
	    if(! country %in% names && nchar(country) <= 3) { # ISO charcode
	        country.code <- match.iso(country)
	        country.idx <- (1:l)[codes==country.code]
	        country.name <- as.character(names[country.idx])
	    } else { # country name
		    country.name <- country
		    country.idx <- (1:l)[names==country.name]
		    country.code <- codes[country.idx]
	    }
		if (length(country.idx) == 0) found <- FALSE
	}
	if (!found) 
		country.name <- country.idx <- country.code <- NULL
	return(list(name=country.name, index=country.idx, code=country.code))
}

match.iso <- local({
    e <- new.env()
    data('iso3166', envir=e)
    get.code <- function(iso) {
        col <- if(nchar(iso) > 2) "charcode3" else "charcode"
        return(e$iso3166$uncode[e$iso3166[[col]] == iso])
    }
})

country.names <- function(meta, countries=NULL, index=FALSE) {
	meta <- get.mcmc.meta(meta)
	if (is.null(countries)) return(meta$regions$country_name)
	if (index) return(meta$regions$country_name[countries])
	get.country.index <- function(code) {
		return(get.country.object(code, meta)$index)
		}
	return(meta$regions$country_name[sapply(countries, get.country.index)])
}

summary.bayesTFR.mcmc <- function(object, country = NULL, par.names = NULL, par.names.cs = NULL, 
							        thin = 1, burnin = 0, ...) {
    res <- list()
    class(res) <- "summary.bayesTFR.mcmc"
    if(is.null(country) && missing(par.names.cs)) par.names.cs <- NULL
	if(is.null(country) && is.null(par.names))
	 	par.names <- if(is.null(object$meta$phase) || object$meta$phase == 2) tfr.parameter.names(trans=TRUE, meta = object$meta) 
	 					else tfr3.parameter.names()
	if(is.null(par.names.cs) && !is.null(country))
		par.names.cs <- if(is.null(object$meta$phase) || object$meta$phase == 2) tfr.parameter.names.cs(trans=TRUE) 
	 					else tfr3.parameter.names.cs()
	if (!is.null(country)) {
		country.obj <- get.country.object(country, object$meta)
		if(is.null(country.obj$name)) stop("Country ", country, " not found.")
		res$country.name <- country.obj$name
		if((is.null(object$meta$phase) || object$meta$phase == 2) && !is.element(country.obj$index, object$meta$id_DL))
			return(res)
		if((!is.null(object$meta$phase) && object$meta$phase == 3) && !is.element(country.obj$index, object$meta$id_phase3))
		    return(res)
		country <- country.obj$code
	} 
	res$results <- summary(coda.mcmc(object, country=country, par.names=par.names,
							par.names.cs=par.names.cs, thin=thin, burnin=burnin), ...)
	return(res)
}

summary.bayesTFR.mcmc.set <- function(object, country=NULL, chain.id=NULL, 
								par.names=NULL, par.names.cs=NULL, 
								meta.only=FALSE, thin=1, burnin=0, ...) {
    if(is.null(country) && missing(par.names.cs)) par.names.cs <- NULL
    res <- list(phase2 = NULL, phase3 = NULL)
	if(is.null(object$meta$phase) || object$meta$phase == 2) {
	  all.standard.names <- c(tfr.parameter.names(meta=object$meta), get.trans.parameter.names(), 
	                          tfr.parameter.names.cs(), get.trans.parameter.names(cs=TRUE), 'tfr')
	  par.names3 <- par.names[par.names %in% tfr3.parameter.names()]
	  par.names3.cs <- par.names.cs[par.names.cs %in% tfr3.parameter.names.cs()]
	  if(is.null(country) && is.null(par.names)) 
		{
		  par.names <- tfr.parameter.names(trans=TRUE, meta = object$meta)
		  if (!is.null(object$mcmc.list[[1]]$uncertainty) && object$mcmc.list[[1]]$uncertainty)
		    par.names3 <- tfr3.parameter.names()
		}
		if(!is.null(country) && is.null(par.names.cs)) 
		{
		  par.names.cs <- tfr.parameter.names.cs(trans=TRUE)
		  if (!is.null(object$mcmc.list[[1]]$uncertainty) && object$mcmc.list[[1]]$uncertainty)
		    par.names3.cs <- tfr3.parameter.names.cs()
		}
	  par.names <- par.names[!(par.names %in% tfr3.parameter.names())]
	  par.names.cs <- par.names.cs[!(par.names.cs %in% tfr3.parameter.names.cs())]
	  
	  if (length(par.names) > 0 || length(par.names.cs) > 0)
	    res$phase2 <- .summary.mcmc.set.phaseII(object, country, chain.id, par.names, par.names.cs, meta.only, thin, burnin, ...)

	  if (length(par.names3) > 0 || length(par.names3.cs) > 0)
	  {
	    if (!is.null(object$mcmc.list[[1]]$uncertainty) && object$mcmc.list[[1]]$uncertainty && has.tfr3.mcmc(object$meta$output.dir))
	    {
	      object3 <- get.tfr3.mcmc(object$meta$output.dir)
	      res$phase3 <- .summary.mcmc.set.phaseIII(object3, country, chain.id, par.names3, par.names3.cs, meta.only, thin, burnin, ...)
	      if (!is.null(res$phase2) && !meta.only && (is.null(country) || (!is.null(country) && res$phase3$estimation)))
	      {
	        stat <- res$phase3$results$statistics
	        quant <- res$phase3$results$quantiles
	        if(is.null(dim(stat))) { # need to add name of the par
	            dim(stat) <- c(1, length(stat))
	            dimnames(stat) <- list(c(par.names3, par.names3.cs), names(res$phase3$results$statistics))
	            dim(quant) <- c(1, length(quant))
	            dimnames(quant) <- list(c(par.names3, par.names3.cs), names(res$phase3$results$quantiles))
	        }
	        if(is.null(dim(res$phase2$results$statistics))) { # need to add name of the par
	          dim(res$phase2$results$statistics) <- c(1, length(res$phase2$results$statistics))
	          dimnames(res$phase2$results$statistics) <- list(c(par.names, par.names.cs), names(res$phase2$results$statistics))
	          dim(res$phase2$results$quantiles) <- c(1, length(res$phase2$results$quantiles))
	          dimnames(res$phase2$results$quantiles) <- list(c(par.names, par.names.cs), names(res$phase2$results$quantiles))
	        }
	        res$phase2$results$statistics <- rbind(res$phase2$results$statistics, stat)
	        res$phase2$results$quantiles <- rbind(res$phase2$results$quantiles, quant)
	        res$phase3$results <- NULL
	      }
	    }
	  }
	} else { # phase III
		if(is.null(country) && is.null(par.names)) par.names <- tfr3.parameter.names()
		if(!is.null(country) && is.null(par.names.cs)) par.names.cs <- tfr3.parameter.names.cs()
		res$phase3 <- .summary.mcmc.set.phaseIII(object, country, chain.id, par.names, par.names.cs, meta.only, thin, burnin, ...)
	}
    class(res) <- "summary.bayesTFR.mcmc.set"
    res
}

.summary.mcmc.set.phaseII <- function(object, country=NULL, chain.id=NULL, 
								par.names=tfr.parameter.names(trans=TRUE), 
								par.names.cs=tfr.parameter.names.cs(trans=TRUE), 
								meta.only=FALSE, thin=1, burnin=0, ...) {
    res <- list(meta = summary(object$meta))
	if(meta.only) 
		return(c(res, list(chain.info = chain.info(object))))
	
	if (!is.null(chain.id))
		return(c(res, list(mcmc = summary(object$mcmc.list[[chain.id]], country=country, par.names=par.names,
							par.names.cs=par.names.cs, thin=thin, burnin=burnin, ...))))
	if (!is.null(country)) {
		country.obj <- get.country.object(country, object$meta)
		if(is.null(country.obj$name)) stop("Country ", country, " not found.")
		res$country.name <- country.obj$name
		if (!is.element(country.obj$index, object$meta$id_DL)) return(res) # not used for estimation
		country <- country.obj$code
	}
	res$results <- summary(coda.list.mcmc(object, country=country, par.names=par.names,
							par.names.cs=par.names.cs, thin=thin, burnin=burnin), ...)
	return(res)
}

.summary.mcmc.set.phaseIII <- function(object, country=NULL, chain.id=NULL, 
								par.names=NULL, par.names.cs=NULL, 
								meta.only=FALSE, thin=1, burnin=0, ...) {
    res <- list(meta = summary(object$meta),
                nr.est = sum(sapply(object$mcmc.list[[1]]$observations, function(x) length(x)-1))
                )
    if(meta.only) return(c(res, list(chain.info = chain.info(object))))
	if (!is.null(chain.id))
	    return(c(res, list(mcmc = summary(object$mcmc.list[[chain.id]], country=country, par.names=par.names,
							par.names.cs=par.names.cs, thin=thin, burnin=burnin, ...))))
	if (!is.null(country)) {
		country.obj <- get.country.object(country, object$meta)
		if(is.null(country.obj$name)) stop("Country ", country, " not found.")
		res$country.name <- country.obj$name
		res$estimation <- FALSE
		if (!is.element(country.obj$index, object$meta$id_phase3)) return(res)
		country <- country.obj$code
		res$estimation <- TRUE
	}
    res$results <- summary(coda.list.mcmc(object, country=country, par.names=par.names,
							par.names.cs=par.names.cs, thin=thin, burnin=burnin), ...)
    return(res)
}

summary.bayesTFR.mcmc.meta <- function(object, ...) {
    res <- list(phase = object$phase, 
                est.period = paste(object$start.year, object$present.year, sep = '-')
                )
    if(is.null(res$phase)) res$phase <- 2
    if(object$phase == 2) {
        res <- c(res, list(nr.countries = object$nr_countries,
                           nr.countries.dl = length(object$id_DL),
                           nr.countries.est = length(object$id_DL[object$id_DL <= object$nr_countries_estimation]),
                           wpp.year = object$wpp.year,
                           annual = if(!is.null(object$annual.simulation)) object$annual.simulation else FALSE
                        ))
    } else { # phase 3
        res <- c(res, list(nr.countries = object$nr.countries, wpp.year = object$parent$wpp.year,
                           annual = if(!is.null(object$parent$annual.simulation)) object$parent$annual.simulation else FALSE))
    }
    class(res) <- "summary.bayesTFR.mcmc.meta"
    return(res)
}

chain.info <- function(object) {
	get.iter <- function(x) x$finished.iter
	res <- list(nr.chains = object$meta$nr.chains, 
				iters = sapply(object$mcmc.list, get.iter),
				thin = object$mcmc.list[[1]]$thin)
	class(res) <- 'summary.chain.info'
	return(res)
}

print.summary.bayesTFR.mcmc.meta <- function(x, ...) {
    if(x$phase == 2) {
        cat('\nMCMCs of phase II')
        cat('\n=================')
        cat('\nNumber of countries:', x$nr.countries)
        if(x$nr.countries.dl < x$nr.countries) 
            cat(" (from which", x$nr.countries - x$nr.countries.dl, "are in phase I)")
        cat('\nHyperparameters estimated using', x$nr.countries.est, 'countries.')
        cat('\nWPP:', x$wpp.year)
        cat('\nInput data: TFR for period', x$est.period)
    } else { # Phase 3
        cat('\nMCMCs of phase III')
        cat('\n==================')
        cat('\nNumber of countries:', x$nr.countries)
        cat('\nWPP:', x$wpp.year)
    }
    cat('\nTime interval: ')
    if(x$annual) cat('annual') else cat('5-year') 
    cat('\n')
}

print.summary.bayesTFR.mcmc.set <- function(x, ...) {
    results <- NULL
    if(!is.null(x$phase2)) {
        p <- x$phase2
        print(p$meta)
        if(!is.null(p$chain.info)) print(p$chain.info)
        if(!is.null(p$mcmc)) print(p$mcmc)
        if(!is.null(p$country.name)){
            cat('\nCountry:', p$country.name, '\n')
            if (is.null(p$results))
                cat('\tnot used for estimation because no decline observed.\n')
        }
        results <- p$results
    }
    if(!is.null(x$phase3)) {
        p <- x$phase3
        print(p$meta)
        cat('Number of observations:', p$nr.est)
        cat('\n')
        if(!is.null(p$chain.info)) print(p$chain.info)
        if(!is.null(p$mcmc)) print(p$mcmc)
        if(!is.null(p$country.name)){
            cat('\nCountry:', p$country.name, '\n')
            if (is.null(p$results) && (!is.null(p$estimation) && !p$estimation || is.null(p$estimation)))
                cat('\tnot used for estimation because it has not reached phase III yet.\n')
        }
        if(!is.null(p$results))
            results <- p$results
    }
    if(!is.null(results)) {
        if(!is.null(x$phase2) && !is.null(x$phase3)) {
            cat('\nParameter estimates for both phases')
            cat('\n===================================')
        }
        print(results)
    }
}

print.bayesTFR.mcmc <- function(x, ...) {
    print(summary(x, ...))
}

print.summary.bayesTFR.mcmc <- function(x, ...) {
    if(!is.null(x$country.name)){
        cat('\nCountry:', x$country.name, '\n')
        if (is.null(x$results))
            cat('\tnot used for estimation.\n')
    }
    if(!is.null(x$results))
        print(x$results)
}

print.bayesTFR.mcmc.set <- function(x, ...) {
    print(summary(x, ...))
}

print.bayesTFR.mcmc.meta <- function(x, ...) {
    print(summary(x, ...))
}

print.bayesTFR.prediction <- function(x, ...) {
    print(summary(x, ...))
}

print.summary.chain.info <- function(x, ...) {
	cat('\nNumber of chains =', x$nr.chains)
	cat('\nIterations =', paste(1,':',sum(x$iters)))
	cat('\nThinning interval =', x$thin)
	cat('\nChains sample sizes:', paste(x$iters, collapse=', '))
	cat('\n')
	invisible(x)
}

get.prediction.summary.data <- function(object, unchanged.pars, country, compact) {
	res <- list()
	for (par in unchanged.pars)
		res[[par]] <- object[[par]]
	proj.and.present.years <- if(is.null(object$proj.years)) 
	      as.integer(dimnames(object$quantiles)[[3]]) else object$proj.years
	res$projection.years <- proj.and.present.years[2:length(proj.and.present.years)]
	if(is.null(country)) return(res)
	if (!is.list(country))
		country <- get.country.object(country, object$mcmc.set$meta)
	if(is.null(country$code)) stop('No prediction available for this country.')
	res$country.name <- country$name
	
	if(compact) {
		quantiles.to.show <- c(0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975)
		quant.index <- is.element(dimnames(object$quantiles)[[2]], as.character(quantiles.to.show))
	} else {
		quant.index <- 1:dim(object$quantiles)[2]
	}
	res$projections <- cbind(t(object$traj.mean.sd[country$index,,]), t(object$quantiles[country$index,quant.index,]))
	colnames(res$projections) <- c('mean', 'SD', 
			paste(as.numeric(dimnames(object$quantiles)[[2]][quant.index])*100, '%', sep=''))
	rownames(res$projections) <- proj.and.present.years
	return(res)
}

.update.summary.data.by.shift <- function(res, object, country) {
	if(is.null(res$projections)) return(res)
	country <- get.country.object(country, object$mcmc.set$meta)
	shift <- get.tfr.shift(country$code, object)
	res$manual <- FALSE
	if(!is.null(shift)) {
		res$projections[,c(1,3:dim(res$projections)[2])] <- res$projections[,c(1,3:dim(res$projections)[2])] + t(matrix(shift,
												 ncol=nrow(res$projections), nrow=ncol(res$projections)-1, byrow=TRUE))
		res$manual <- TRUE
	}
	return(res)
}

summary.bayesTFR.prediction <- function(object, country=NULL, compact=TRUE, ...) {
	res <- get.prediction.summary.data(object, 
				unchanged.pars=c('burnin', 'thin', 'nr.traj', 'mu', 'rho', 'sigmaAR1', 'use.tfr3', 'burnin3', 'thin3'), 
				country=country, compact=compact)
	class(res) <- 'summary.bayesTFR.prediction'
	return(.update.summary.data.by.shift(res, object, country))
}

print.summary.bayesTFR.prediction <- function(x, digits = 3, ...) {
	cat('\nProjections:', length(x$projection.years), '(', x$projection.years[1], '-', 
					x$projection.years[length(x$projection.years)], ')')
	cat('\nTrajectories:', x$nr.traj)
	cat('\nPhase II burnin:', x$burnin)
	cat('\nPhase II thin:', x$thin)
	if(x$use.tfr3) {
		cat('\nPhase III burnin:', x$burnin3)
		cat('\nPhase III thin:', x$thin3)
		cat('\n')
	} else {
		cat('\nParameters of AR(1):\n')
		arpars <- c(x$mu, x$rho, x$sigmaAR1)
		ar.table<-matrix(arpars, nrow=1, ncol=length(arpars),
						dimnames=list('',c('mu', 'rho', rep('sigma', length(x$sigmaAR1)))))
		print(ar.table, digits=digits, ...)
	}
	if(!is.null(x$country.name)) {
		cat('\nCountry:', x$country.name, '\n')
		cat('\nProjected TFR:')
		if(x$manual) cat(' (values have been manually modified)')
		cat('\n')
		print(x$projections, digits=digits, ...)
	}
}


summary.bayesTFR.convergence <- function(object, expand = FALSE, ...) {
    res <- list(burnin = object$burnin, thin = object$thin, express = object$express,
                results = summary(object$mcmc.set, meta.only=TRUE))
	if(!is.null(object$thin.mcmc)) 
	    res$thin.mcmc <- chain.info(object$thin.mcmc)
	if(!object$express) {
	    res$nr.countries.used <- object$nr.countries['used']
	    res$nr.countries.total <- object$nr.countries['total']
	}
	if(nrow(object$result) > 0) {
		if (object$lresult.country.independent > 0) 
		    res$hyper.pars <- object$result[1:object$lresult.country.independent,,drop=FALSE]
		if (nrow(object$result) - object$lresult.country.independent > 0)
		    res$cntry.pars <- object$result[(object$lresult.country.independent+1):nrow(object$result),,drop=FALSE]
		res$iter.needed <- object$iter.needed
	}
	res$status <- names(object$status)[object$status]
	if(object$status['green']) 
	    res$use.nr.traj <- object$use.nr.traj

	if(expand && !object$status['green'] && !is.null(object$country.specific) 
			&& !is.null(object$country.specific$not.converged.parameters)) 
	    res$not.converged.pars <- object$country.specific$not.converged.parameters
	class(res) <- "summary.bayesTFR.convergence"
	return(res)
}

print.summary.bayesTFR.convergence <- function(x, ...) {
    cat('\nConvergence diagnostics for burnin =', x$burnin, 'and thin =', x$thin)
    cat('\n********************************************************')
    cat('\nFull chains:')
    cat('\n============')
    print(x$results)
    
    if(!is.null(x$thin.mcmc)) {
        cat('\nThinned, burned and collapsed chain:')
        cat('\n====================================')
        print(x$thin.mcmc)
    }
    if(x$express) cat('\nExpress diagnostics - no country-specific parameters were included.')
    else cat('\nConvergence checked on', x$nr.countries.used, 
             'countries out of', x$nr.countries.total, 'countries total.')
    cat('\n')
    if(!is.null(x$iter.needed)) {
        if (!is.null(x$hyper.pars)) {
            cat('\nNot converged country-independent parameters:')
            cat('\n---------------------------------------------\n')
            print(x$hyper.pars)
        }
        if (!is.null(x$cntry.pars)) {
            cat('\nNot converged country-specific parameters:')
            cat('\n-------------------------------------------\n')
            print(x$cntry.pars)
        }
        cat('\n\nAt least', x$iter.needed, 'more iterations needed  to achieve convergence.\n')
    }
    if(x$status == "green") {
        cat('\nSimulation has converged.')
        cat('\nNumber of trajectories to be used:', x$use.nr.traj)
    }
    if(!is.null(x$not.converged.pars)) {
        cat('\nWarning: The following parameters did not converge:\n')
        print(x$not.converged.pars)
    }
    cat('\nStatus:', x$status, '\n')
}

print.bayesTFR.convergence <- function(x, ...) 
    print(summary(x))

tfr.info <- function(sim.dir) {
	mc <- get.tfr.mcmc(sim.dir=sim.dir)
	if (is.null(mc)) {
		cat('No TFR simulation available in', sim.dir)
		return (NULL)
	}
	summary(mc)
}

get.cov.gammas <- function(mcmc.set=NULL, sim.dir = NULL, burnin = 200, chain.id=1){
	# this is for one chain
	if (is.null(mcmc.set))
		mcmc.set <- get.tfr.mcmc(sim.dir=sim.dir)
 	cov_gammas_cii = array(NA, c(mcmc.set$meta$nr_countries, 3,3))
 	for (country in mcmc.set$meta$id_DL){
 		country.obj <- get.country.object(country, mcmc.set$meta, index=TRUE)
  		gamma_si <- load.tfr.parameter.traces.cs(mcmc.set$mcmc.list[[chain.id]], 
  							par.names='gamma', country=country.obj$code, burnin=burnin)
        cov_gammas_cii[country,, ] <- cov(gamma_si)
 	}
 	cov_gammas_cii = (2.4^2/3)*cov_gammas_cii 
 	return(list(values=cov_gammas_cii, country_codes=mcmc.set$meta$regions$country_code))
}

get.thinning.index <- function(nr.points, all.points) {
 	if (!is.null(nr.points)) {
		nr.points <- ifelse(nr.points >= all.points, all.points, nr.points)
	} else {
		nr.points <- all.points
	}
	if (nr.points > 0) {
		step <- all.points/nr.points
		idx <- floor(seq(floor(step), all.points, by=step))
	} else idx<-NULL
	return(list(nr.points=nr.points, index=idx))
}

burn.and.thin <- function(mcmc, burnin=0, thin=1) {
	# Return thin and burnin that is consolidated with the original thin usd for storing mcmc.
	# If there is need for more thinning, it returns the corresponding thinning index
	th.burnin <- get.thinned.burnin(mcmc,burnin)
	index <- NULL
	if (thin > mcmc$thin) {
    	index <- unique(round(seq(1,mcmc$length-th.burnin, by=thin/mcmc$thin)))
    }
    thin <- max(mcmc$thin, thin) # thin cannot be smaller than the original thin
	return(list(thin=thin, burnin=th.burnin, index=index))
}

no.traces.loaded <- function(mcmc) return((length(mcmc$traces) == 1) && mcmc$traces == 0)

tfr.set.identical <- function(mcmc.set1, mcmc.set2, include.output.dir=TRUE) {
	# Test if two bayesTFR sets are identical
	if(!include.output.dir) 
		mcmc.set1$meta$output.dir <- mcmc.set2$meta$output.dir <- NULL

	same <- setequal(names(mcmc.set1), names(mcmc.set2)) && tfr.meta.identical(mcmc.set1$meta, mcmc.set2$meta) && length(mcmc.set1$mcmc.list) == length(mcmc.set2$mcmc.list)
	if(!same) return(same)
	for(i in 1:length(mcmc.set1$mcmc.list)) {
		if(!include.output.dir) mcmc.set1$mcmc.list[[i]]$meta$output.dir <- mcmc.set2$mcmc.list[[i]]$meta$output.dir <- NULL
		same <- same && tfr.identical(mcmc.set1$mcmc.list[[i]], mcmc.set2$mcmc.list[[i]])
	}
	return(same)
}

tfr.meta.identical <- function(meta1, meta2) {
    identical(meta1, meta2) || all.equal(meta1, meta2)
}

tfr.identical <- function(mcmc1, mcmc2) {
	# Test if two mcmcs are identical
	same <- setequal(names(mcmc1), names(mcmc2))
	if(!same || !tfr.meta.identical(mcmc1$meta, mcmc2$meta)) return(FALSE)
	for(item in setdiff(names(mcmc1), "meta")) 
		same <- same && identical(mcmc1[[item]], mcmc2[[item]])
	return(same)
}

get.tfr.trajectories <- function(tfr.pred, country) {
	country.obj <- get.country.object(country, tfr.pred$mcmc.set$meta)
	return(get.trajectories(tfr.pred, country.obj$code)$trajectories)
}

"get.nr.countries" <- function(meta, ...) UseMethod("get.nr.countries")
 
get.nr.countries.bayesTFR.mcmc.meta <- function(meta, ...) 
	return (if(is.null(meta$phase) || (meta$phase==2)) meta$nr_countries else meta$nr.countries)

"get.nrest.countries" <- function(meta, ...) UseMethod("get.nrest.countries")
 
get.nrest.countries.bayesTFR.mcmc.meta <- function(meta, ...) 
	return (if(is.null(meta$phase) || (meta$phase==2)) meta$nr_countries_estimation else meta$nr.countries)

"get.data.matrix" <- function(meta, ...) UseMethod("get.data.matrix")
 
get.data.matrix.bayesTFR.mcmc.meta <- function(meta, ...) return (meta$tfr_matrix)

"get.countries.index" <- function(meta, ...) UseMethod("get.countries.index")

get.countries.index.bayesTFR.mcmc.meta  <- function(meta, ...) 
	return (if(is.null(meta$phase) || (meta$phase==2)) meta$id_DL else meta$id_phase3)

"get.countries.table" <- function(object, ...) UseMethod("get.countries.table")
get.countries.table.bayesTFR.mcmc.set <- function(object, iso = FALSE, ...) {
	ctable <- data.frame(code=object$meta$regions$country_code, name=object$meta$regions$country_name)
	if(!is.null(object$meta$phase) && (object$meta$phase==3)) ctable <- ctable[object$meta$id_phase3,]
	if(iso) ctable <- add.iso(ctable)
	return(ctable)
}

get.countries.table.bayesTFR.prediction <- function(object, iso = FALSE, ...) {
	n <- dim(get.data.imputed(object))[2]
	ctable <- data.frame(code=object$mcmc.set$meta$regions$country_code[1:n], 
					name=object$mcmc.set$meta$regions$country_name[1:n])
	if(iso) ctable <- add.iso(ctable)
	return(ctable)
}

add.iso <- function(tbl){
    iso <- get(data('iso3166', envir=environment()))[, c("uncode", "charcode", "charcode3")]
    colnames(iso) <- c("code", "iso2", "iso3")
    merge(tbl, iso, by = "code", sort = FALSE, all.x = TRUE)
}

"get.countries.phase" <- function(object, ...) UseMethod("get.countries.phase")
get.countries.phase.bayesTFR.mcmc.set <- function(object, ...) {
    ctable <- data.frame(code=object$meta$regions$country_code, name=object$meta$regions$country_name,
                         phase = 1)
    ctable$phase[object$meta$id_DL] <- 2 # ever experienced Phase II
    ctable$phase[which(object$meta$lambda_c < object$meta$T_end_c)] <- 3 # now in phase III
    return(ctable)
}

get.countries.phase.bayesTFR.prediction <- function(object, ...) 
    get.countries.phase(object$mcmc.set)

get.item <- function(l, item, default)
    ifelse(item %in% names(l), l[[item]], default)