diff --git a/R/bridge_sampler.R b/R/bridge_sampler.R index 0c1335f..f7e2c54 100644 --- a/R/bridge_sampler.R +++ b/R/bridge_sampler.R @@ -263,9 +263,10 @@ bridge_sampler.stanfit <- function(samples = NULL, stanfit_model = samples, param_types = rep("real", ncol(samples_4_fit)), transTypes = transTypes, repetitions = repetitions, cores = cores, - packages = "rstan", maxiter = maxiter, silent = silent, - verbose = verbose, - r0 = 0.5, tol1 = 1e-10, tol2 = 1e-4)) + packages = c("rstan", "mvtnorm"), + maxiter = maxiter, silent = silent, + verbose = verbose, r0 = 0.5, tol1 = 1e-10, + tol2 = 1e-4)) } else { bridge_output <- do.call(what = paste0(".bridge.sampler.", method), args = list(samples_4_fit = samples_4_fit, @@ -278,9 +279,10 @@ bridge_sampler.stanfit <- function(samples = NULL, stanfit_model = samples, transTypes = transTypes, repetitions = repetitions, varlist = "stanfit", envir = sys.frame(sys.nframe()), - cores = cores, packages = "rstan", maxiter = maxiter, - silent = silent, verbose = verbose, - r0 = 0.5, tol1 = 1e-10, tol2 = 1e-4)) + cores = cores, packages = c("rstan", "mvtnorm"), + maxiter = maxiter, silent = silent, + verbose = verbose, r0 = 0.5, tol1 = 1e-10, + tol2 = 1e-4)) } return(bridge_output) @@ -490,7 +492,7 @@ bridge_sampler.stanreg <- bridge_output <- bridge_sampler(samples = samples, log_posterior = .stan_log_posterior, data = list(stanfit = sf), lb = lb, ub = ub, repetitions = repetitions, method = method, cores = cores, - use_neff = use_neff, packages = "rstan", + use_neff = use_neff, packages = c("rstan", "mvtnorm"), maxiter = maxiter, silent = silent, verbose = verbose) } else { @@ -500,7 +502,7 @@ bridge_sampler.stanreg <- repetitions = repetitions, varlist = "stanfit", envir = sys.frame(sys.nframe()), method = method, cores = cores, use_neff = use_neff, - packages = "rstan", maxiter = maxiter, + packages = c("rstan", "mvtnorm"), maxiter = maxiter, silent = silent, verbose = verbose) } return(bridge_output) diff --git a/R/bridge_sampler_normal.R b/R/bridge_sampler_normal.R index 6cf550c..72873d8 100644 --- a/R/bridge_sampler_normal.R +++ b/R/bridge_sampler_normal.R @@ -38,12 +38,15 @@ # sample from multivariate normal distribution and evaluate for posterior samples and generated samples q12 <- dmvnorm(samples_4_iter, mean = m, sigma = V, log = TRUE) - gen_samples <- vector(mode = "list", length = repetitions) - q22 <- vector(mode = "list", length = repetitions) - for (i in seq_len(repetitions)) { - gen_samples[[i]] <- rmvnorm(n_post, mean = m, sigma = V) - colnames(gen_samples[[i]]) <- colnames(samples_4_iter) - q22[[i]] <- dmvnorm(gen_samples[[i]], mean = m, sigma = V, log = TRUE) + if (cores == 1 || repetitions == 1) { + # if cores > 1 and repetitions > 1 we do this below using parallelization + gen_samples <- vector(mode = "list", length = repetitions) + q22 <- vector(mode = "list", length = repetitions) + for (i in seq_len(repetitions)) { + gen_samples[[i]] <- rmvnorm(n_post, mean = m, sigma = V) + colnames(gen_samples[[i]]) <- colnames(samples_4_iter) + q22[[i]] <- dmvnorm(gen_samples[[i]], mean = m, sigma = V, log = TRUE) + } } # evaluate log of likelihood times prior for posterior samples and generated samples @@ -57,6 +60,19 @@ } } else if (cores > 1) { if ( .Platform$OS.type == "unix") { + # sample from multivariate normal distribution and evaluate for posterior samples and generated samples + if (repetitions > 1) { + gen_samples <- parallel::mclapply(seq_len(repetitions), FUN = + function(x) rmvnorm(n_post, mean = m, sigma = V), + mc.preschedule = FALSE, + mc.cores = cores) + sapply(seq_along(gen_samples), function(i) colnames(gen_samples[[i]]) <- colnames(samples_4_iter)) + q22 <- parallel::mclapply(seq_along(gen_samples), FUN = + function(i) dmvnorm(gen_samples[[i]], mean = m, sigma = V, log = TRUE), + mc.preschedule = FALSE, + mc.cores = cores) + } + split1 <- .split_matrix(matrix=.invTransform2Real(samples_4_iter, lb, ub, param_types), cores=cores) q11 <- parallel::mclapply(split1, FUN = function(x) apply(x, 1, log_posterior, data = data, ...), @@ -73,8 +89,23 @@ } } else { cl <- parallel::makeCluster(cores, useXDR = FALSE) + on.exit(parallel::stopCluster(cl)) sapply(packages, function(x) parallel::clusterCall(cl = cl, "require", package = x, character.only = TRUE)) + + if (repetitions > 1) { + # sample from multivariate normal distribution and evaluate for posterior samples and generated samples + parallel::clusterExport(cl = cl, varlist = c("m", "V", "n_post"), envir = environment()) + gen_samples <- parallel::clusterApplyLB(cl = cl, + x = seq_len(repetitions), + fun = function(i) rmvnorm(n_post, mean = m, sigma = V)) + sapply(seq_along(gen_samples), function(i) colnames(gen_samples[[i]]) <- colnames(samples_4_iter)) + parallel::clusterExport(cl = cl, varlist = "gen_samples", envir = environment()) + q22 <- parallel::clusterApplyLB(cl = cl, + x = seq_along(gen_samples), + fun = function(i) dmvnorm(gen_samples[[i]], mean = m, sigma = V, log = TRUE)) + } + parallel::clusterExport(cl = cl, varlist = varlist, envir = envir) if ( ! is.null(rcppFile)) { @@ -91,7 +122,7 @@ q21[[i]] <- parallel::parRapply(cl = cl, x = .invTransform2Real(gen_samples[[i]], lb, ub, param_types), log_posterior, data = data, ...) + .logJacobian(gen_samples[[i]], transTypes, lb, ub) } - parallel::stopCluster(cl) + # parallel::stopCluster(cl) } } if(verbose) { diff --git a/R/bridge_sampler_warp3.R b/R/bridge_sampler_warp3.R index 559e51a..5b20587 100644 --- a/R/bridge_sampler_warp3.R +++ b/R/bridge_sampler_warp3.R @@ -39,21 +39,23 @@ # sample from multivariate normal distribution and evaluate for posterior samples and generated samples q12 <- dmvnorm((samples_4_iter - matrix(m, nrow = n_post, ncol = length(m), byrow = TRUE)) %*% - t(solve(L)), sigma = diag(ncol(samples_4_fit)), log = TRUE) - q22 <- vector(mode = "list", length = repetitions) - gen_samples <- vector(mode = "list", length = repetitions) - for (i in seq_len(repetitions)) { - gen_samples[[i]] <- rmvnorm(n_post, sigma = diag(ncol(samples_4_fit))) - colnames(gen_samples[[i]]) <- colnames(samples_4_iter) - q22[[i]] <- dmvnorm(gen_samples[[i]], sigma = diag(ncol(samples_4_fit)), log = TRUE) + t(solve(L)), sigma = diag(ncol(samples_4_fit)), log = TRUE) + if (cores == 1 || repetitions == 1) { + # if cores > 1 or repetitions > 1, we do this below using parallelization + q22 <- vector(mode = "list", length = repetitions) + gen_samples <- vector(mode = "list", length = repetitions) + for (i in seq_len(repetitions)) { + gen_samples[[i]] <- rmvnorm(n_post, sigma = diag(ncol(samples_4_fit))) + colnames(gen_samples[[i]]) <- colnames(samples_4_iter) + q22[[i]] <- dmvnorm(gen_samples[[i]], sigma = diag(ncol(samples_4_fit)), log = TRUE) + } } e <- as.brob( exp(1) ) - + # evaluate log of likelihood times prior for posterior samples and generated samples q21 <- vector(mode = "list", length = repetitions) if (cores == 1) { - q11 <- log(e^(apply(.invTransform2Real(samples_4_iter, lb, ub, param_types), 1, log_posterior, data = data,...) + .logJacobian(samples_4_iter, transTypes, lb, ub)) + e^(apply(.invTransform2Real(matrix(2*m, nrow = n_post, ncol = length(m), byrow = TRUE) - @@ -72,6 +74,18 @@ } } else if (cores > 1) { if ( .Platform$OS.type == "unix") { + if (repetitions > 1) { + gen_samples <- parallel::mclapply(seq_len(repetitions), FUN = + function(x) rmvnorm(n_post, sigma = diag(ncol(samples_4_fit))), + mc.preschedule = FALSE, + mc.cores = cores) + sapply(seq_along(gen_samples), function(i) colnames(gen_samples[[i]]) <- colnames(samples_4_iter)) + q22 <- parallel::mclapply(seq_along(gen_samples), FUN = + function(i) dmvnorm(gen_samples[[i]], sigma = diag(ncol(samples_4_fit)), log = TRUE), + mc.preschedule = FALSE, + mc.cores = cores) + } + split1a <- .split_matrix(matrix=.invTransform2Real(samples_4_iter, lb, ub, param_types), cores=cores) split1b <- .split_matrix(matrix=.invTransform2Real( matrix(2*m, nrow = n_post, ncol = length(m), byrow = TRUE) - samples_4_iter, lb, ub, param_types @@ -113,9 +127,23 @@ tmp_mat2, transTypes, lb, ub))) } } else { + # sample from multivariate normal distribution and evaluate for posterior samples and generated samples cl <- parallel::makeCluster(cores, useXDR = FALSE) + on.exit(parallel::stopCluster(cl)) sapply(packages, function(x) parallel::clusterCall(cl = cl, "require", package = x, character.only = TRUE)) - + + if (repetitions > 1) { + parallel::clusterExport(cl = cl, varlist = c("n_post", "samples_4_fit"), envir = environment()) + gen_samples <- parallel::clusterApplyLB(cl = cl, + x = seq_len(repetitions), + fun = function(x) rmvnorm(n_post, sigma = diag(ncol(samples_4_fit)))) + sapply(seq_along(gen_samples), function(i) colnames(gen_samples[[i]]) <- colnames(samples_4_iter)) + parallel::clusterExport(cl = cl, varlist = "gen_samples", envir = environment()) + q22 <- parallel::clusterApplyLB(cl = cl, + x = seq_along(gen_samples), + fun = function(i) dmvnorm(gen_samples[[i]], sigma = diag(ncol(samples_4_fit)), log = TRUE)) + } + parallel::clusterExport(cl = cl, varlist = varlist, envir = envir) if ( ! is.null(rcppFile)) { @@ -152,7 +180,7 @@ .logJacobian(matrix(m, nrow = n_post, ncol = length(m), byrow = TRUE) + gen_samples[[i]] %*% t(L), transTypes, lb, ub))) } - parallel::stopCluster(cl) + # parallel::stopCluster(cl) } }