R/RcppExports.R

In spBPS: Bayesian Predictive Stacking for Scalable Geospatial Transfer Learning

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' Compute the BPS weights by convex optimization
#'
#' @param scores [matrix] \eqn{N \times K} of expected predictive density evaluations for the K models considered
#'
#' @return conv_opt [function] to perform convex optimiazion with CVXR R package
#'
CVXR_opt <- function(scores) {
    .Call(`_spBPS_CVXR_opt`, scores)
}

#' Combine subset models wiht BPS
#'
#' @param fit_list [list] K fitted model outputs composed by two elements each: first named \eqn{epd}, second named \eqn{W}
#' @param K [integer] number of folds
#' @param rp [double] percentage of observations to take into account for optimization (\code{default=1})
#'
#' @return [matrix] posterior predictive density evaluations (each columns represent a different model)
#'
#' @examples
#' \donttest{
#' ## Generate subsets of data
#' n <- 100
#' p <- 3
#' X <- matrix(rnorm(n*p), nrow = n, ncol = p)
#' Y <- matrix(rnorm(n), nrow = n, ncol = 1)
#' crd <- matrix(runif(n*2), nrow = n, ncol = 2)
#' data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10)
#'
#' ## Select competitive set of values for hyperparameters
#' delta_seq <- c(0.1, 0.2, 0.3)
#' phi_seq <- c(3, 4, 5)
#'
#' ## Perform Bayesian Predictive Stacking within subsets
#' fit_list <- vector(length = 10, mode = "list")
#' for (i in 1:10) {
#'     Yi <- data_part$Y_list[[i]]
#'     Xi <- data_part$X_list[[i]]
#'     crd_i <- data_part$crd_list[[i]]
#'     p <- ncol(Xi)
#'     bps <- spBPS::BPS_weights(data = list(Y = Yi, X = Xi),
#'                                priors = list(mu_b = matrix(rep(0, p)),
#'                                              V_b = diag(10, p),
#'                                              a = 2,
#'                                              b = 2), coords = crd_i,
#'                                              hyperpar = list(delta = delta_seq,
#'                                                              phi = phi_seq),
#'                                              K = 5)
#'      w_hat <- bps$W
#'      epd <- bps$epd
#'      fit_list[[i]] <- list(epd, w_hat) }
#'
#' ## Combination weights between partitions using Bayesian Predictive Stacking
#' comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1)
#' }
#'
#' @export
BPS_combine <- function(fit_list, K, rp) {
    .Call(`_spBPS_BPS_combine`, fit_list, K, rp)
}

#' Combine subset models wiht Pseudo-BMA
#'
#' @param fit_list [list] K fitted model outputs composed by two elements each: first named \eqn{epd}, second named \eqn{W}
#'
#' @return [matrix] posterior predictive density evaluations (each columns represent a different model)
#'
#' @examples
#' \donttest{
#' ## Generate subsets of data
#' n <- 100
#' p <- 3
#' X <- matrix(rnorm(n*p), nrow = n, ncol = p)
#' Y <- matrix(rnorm(n), nrow = n, ncol = 1)
#' crd <- matrix(runif(n*2), nrow = n, ncol = 2)
#' data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10)
#'
#' ## Select competitive set of values for hyperparameters
#' delta_seq <- c(0.1, 0.2, 0.3)
#' phi_seq <- c(3, 4, 5)
#'
#' ## Perform Bayesian Predictive Stacking within subsets
#' fit_list <- vector(length = 10, mode = "list")
#' for (i in 1:10) {
#'     Yi <- data_part$Y_list[[i]]
#'     Xi <- data_part$X_list[[i]]
#'     crd_i <- data_part$crd_list[[i]]
#'     p <- ncol(Xi)
#'     bps <- spBPS::BPS_weights(data = list(Y = Yi, X = Xi),
#'                                priors = list(mu_b = matrix(rep(0, p)),
#'                                              V_b = diag(10, p),
#'                                              a = 2,
#'                                              b = 2), coords = crd_i,
#'                                              hyperpar = list(delta = delta_seq,
#'                                                              phi = phi_seq),
#'                                              K = 5)
#'      w_hat <- bps$W
#'      epd <- bps$epd
#'      fit_list[[i]] <- list(epd, w_hat) }
#'
#' ## Combination weights between partitions using Pseudo Bayesian Model Averaging
#' comb_bps <- BPS_PseudoBMA(fit_list = fit_list)
#' }
#'
#' @export
BPS_PseudoBMA <- function(fit_list) {
    .Call(`_spBPS_BPS_PseudoBMA`, fit_list)
}

#' Compute the parameters for the posteriors distribution of \eqn{\beta} and \eqn{\Sigma} (i.e. updated parameters)
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param priors [list] priors: named \eqn{\mu_b},\eqn{V_b},\eqn{a},\eqn{b}
#' @param coords [matrix] sample coordinates for X and Y
#' @param hyperpar [list] two elemets: first named \eqn{\delta}, second named \eqn{\phi}
#'
#' @return [list] posterior update parameters
#'
fit_cpp <- function(data, priors, coords, hyperpar) {
    .Call(`_spBPS_fit_cpp`, data, priors, coords, hyperpar)
}

#' Sample R draws from the posterior distributions
#'
#' @param poster [list] output from \code{fit_cpp} function
#' @param R [integer] number of posterior samples
#' @param par if \code{TRUE} only \eqn{\beta} and \eqn{\sigma^2} are sampled (\eqn{\omega} is omitted)
#' @param p [integer] if \code{par = TRUE}, it specifies the column number of \eqn{X}
#'
#' @return [list] posterior samples
#'
post_draws <- function(poster, R, par, p) {
    .Call(`_spBPS_post_draws`, poster, R, par, p)
}

#' Draw from the joint posterior predictive for a set of unobserved covariates
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param X_u [matrix] unobserved instances covariate matrix
#' @param d_u [matrix] unobserved instances distance matrix
#' @param d_us [matrix] cross-distance between unobserved and observed instances matrix
#' @param hyperpar [list] two elemets: first named \eqn{\delta}, second named \eqn{\phi}
#' @param poster [list] output from \code{fit_cpp} function
#' @param R [integer] number of posterior predictive samples
#'
#' @return [list] posterior predictive samples
#'
r_pred_joint <- function(data, X_u, d_u, d_us, hyperpar, poster, R) {
    .Call(`_spBPS_r_pred_joint`, data, X_u, d_u, d_us, hyperpar, poster, R)
}

#' Draw from the marginals posterior predictive for a set of unobserved covariates
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param X_u [matrix] unobserved instances covariate matrix
#' @param d_u [matrix] unobserved instances distance matrix
#' @param d_us [matrix] cross-distance between unobserved and observed instances matrix
#' @param hyperpar [list] two elemets: first named \eqn{\delta}, second named \eqn{\phi}
#' @param poster [list] output from \code{fit_cpp} function
#' @param R [integer] number of posterior predictive samples
#'
#' @return [list] posterior predictive samples
#'
r_pred_marg <- function(data, X_u, d_u, d_us, hyperpar, poster, R) {
    .Call(`_spBPS_r_pred_marg`, data, X_u, d_u, d_us, hyperpar, poster, R)
}

#' Draw from the conditional posterior predictive for a set of unobserved covariates
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param X_u [matrix] unobserved instances covariate matrix
#' @param d_u [matrix] unobserved instances distance matrix
#' @param d_us [matrix] cross-distance between unobserved and observed instances matrix
#' @param hyperpar [list] two elemets: first named \eqn{\delta}, second named \eqn{\phi}
#' @param poster [list] output from \code{fit_cpp} function
#' @param post [list] output from \code{post_draws} function
#'
#' @return [list] posterior predictive samples
#'
r_pred_cond <- function(data, X_u, d_u, d_us, hyperpar, poster, post) {
    .Call(`_spBPS_r_pred_cond`, data, X_u, d_u, d_us, hyperpar, poster, post)
}

#' Evaluate the density of a set of unobserved response with respect to the conditional posterior predictive
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param X_u [matrix] unobserved instances covariate matrix
#' @param Y_u [matrix] unobserved instances response matrix
#' @param d_u [matrix] unobserved instances distance matrix
#' @param d_us [matrix] cross-distance between unobserved and observed instances matrix
#' @param hyperpar [list] two elemets: first named \eqn{\delta}, second named \eqn{\phi}
#' @param poster [list] output from \code{fit_cpp} function
#'
#' @return [vector] posterior predictive density evaluations
#'
d_pred_cpp <- function(data, X_u, Y_u, d_u, d_us, hyperpar, poster) {
    .Call(`_spBPS_d_pred_cpp`, data, X_u, Y_u, d_u, d_us, hyperpar, poster)
}

#' Compute the LOOCV of the density evaluations for fixed values of the hyperparameters
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param priors [list] priors: named \eqn{\mu_b},\eqn{V_b},\eqn{a},\eqn{b}
#' @param coords [matrix] sample coordinates for X and Y
#' @param hyperpar [list] two elemets: first named \eqn{\delta}, second named \eqn{\phi}
#'
#' @return [vector] posterior predictive density evaluations
#'
dens_loocv <- function(data, priors, coords, hyperpar) {
    .Call(`_spBPS_dens_loocv`, data, priors, coords, hyperpar)
}

#' Compute the KCV of the density evaluations for fixed values of the hyperparameters
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param priors [list] priors: named \eqn{\mu_b},\eqn{V_b},\eqn{a},\eqn{b}
#' @param coords [matrix] sample coordinates for X and Y
#' @param hyperpar [list] two elemets: first named \eqn{\delta}, second named \eqn{\phi}
#' @param K [integer] number of folds
#'
#' @return [vector] posterior predictive density evaluations
#'
dens_kcv <- function(data, priors, coords, hyperpar, K) {
    .Call(`_spBPS_dens_kcv`, data, priors, coords, hyperpar, K)
}

#' Return the CV predictive density evaluations for all the model combinations
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param priors [list] priors: named \eqn{\mu_b},\eqn{V_b},\eqn{a},\eqn{b}
#' @param coords [matrix] sample coordinates for X and Y
#' @param hyperpar [list] two elemets: first named \eqn{\delta}, second named \eqn{\phi}
#' @param useKCV if \code{TRUE} K-fold cross validation is used instead of LOOCV (no \code{default})
#' @param K [integer] number of folds
#'
#' @return [matrix] posterior predictive density evaluations (each columns represent a different model)
#'
models_dens <- function(data, priors, coords, hyperpar, useKCV, K) {
    .Call(`_spBPS_models_dens`, data, priors, coords, hyperpar, useKCV, K)
}

#' Compute the BPS weights by convex optimization
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param priors [list] priors: named \eqn{\mu_b},\eqn{V_b},\eqn{a},\eqn{b}
#' @param coords [matrix] sample coordinates for X and Y
#' @param hyperpar [list] two elemets: first named \eqn{\delta}, second named \eqn{\phi}
#' @param K [integer] number of folds
#'
#' @return [matrix] posterior predictive density evaluations (each columns represent a different model)
#'
#' @examples
#' ## Generate subsets of data
#' n <- 100
#' p <- 3
#' X <- matrix(rnorm(n*p), nrow = n, ncol = p)
#' Y <- matrix(rnorm(n), nrow = n)
#' crd <- matrix(runif(n*2), nrow = n, ncol = 2)
#'
#' ## Select competitive set of values for hyperparameters
#' delta_seq <- c(0.1, 0.2, 0.3)
#' phi_seq <- c(3, 4, 5)
#'
#' ## Perform Bayesian Predictive Stacking within subsets
#' bps <- spBPS::BPS_weights(data = list(Y = Y, X = X),
#'                                priors = list(mu_b = matrix(rep(0, p)),
#'                                              V_b = diag(10, p),
#'                                              a = 2,
#'                                              b = 2), coords = crd,
#'                                              hyperpar = list(delta = delta_seq,
#'                                                              phi = phi_seq),
#'                                              K = 5)
#'
#' @export
BPS_weights <- function(data, priors, coords, hyperpar, K) {
    .Call(`_spBPS_BPS_weights`, data, priors, coords, hyperpar, K)
}

#' Compute the BPS spatial prediction given a set of stacking weights
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param X_u [matrix] unobserved instances covariate matrix
#' @param priors [list] priors: named \eqn{\mu_b},\eqn{V_b},\eqn{a},\eqn{b}
#' @param coords [matrix] sample coordinates for X and Y
#' @param crd_u [matrix] unboserved instances coordinates
#' @param hyperpar [list] two elemets: first named \eqn{\delta}, second named \eqn{\phi}
#' @param W [matrix] set of stacking weights
#' @param R [integer] number of desired samples
#'
#' @return [list] BPS posterior predictive samples
#'
#' @examples
#' \donttest{
#' ## Generate subsets of data
#' n <- 100
#' p <- 3
#' X <- matrix(rnorm(n*p), nrow = n, ncol = p)
#' Y <- matrix(rnorm(n), nrow = n, ncol = 1)
#' crd <- matrix(runif(n*2), nrow = n, ncol = 2)
#' data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10)
#'
#' ## Select competetive set of values for hyperparameters
#' delta_seq <- c(0.1, 0.2, 0.3)
#' phi_seq <- c(3, 4, 5)
#'
#' ## Fit local models
#' fit_list <- vector(length = 10, mode = "list")
#' for (i in 1:10) {
#'     Yi <- data_part$Y_list[[i]]
#'     Xi <- data_part$X_list[[i]]
#'     crd_i <- data_part$crd_list[[i]]
#'     p <- ncol(Xi)
#'     bps <- spBPS::BPS_weights(data = list(Y = Yi, X = Xi),
#'                                priors = list(mu_b = matrix(rep(0, p)),
#'                                              V_b = diag(10, p),
#'                                              a = 2,
#'                                              b = 2), coords = crd_i,
#'                                              hyperpar = list(delta = delta_seq,
#'                                                              phi = phi_seq),
#'                                              K = 5)
#'      w_hat <- bps$W
#'      epd <- bps$epd
#'      fit_list[[i]] <- list(epd, w_hat) }
#'
#' ## Model combination weights between partitions using Bayesian Predictive Stacking
#' comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1)
#' Wbps <- comb_bps$W
#' W_list <- comb_bps$W_list
#'
#' ## Generate prediction points
#' m <- 50
#' X_new <- matrix(rnorm(m*p), nrow = m, ncol = p)
#' crd_new <- matrix(runif(m*2), nrow = m, ncol = 2)
#'
#' ## Perform posterior predictive sampling
#' R <- 250
#' subset_ind <- sample(1:10, R, TRUE, Wbps)
#' predictions <- vector(length = R, mode = "list")
#' for (r in 1:R) {
#'   ind_s <- subset_ind[r]
#'   Ys <- matrix(data_part$Y_list[[ind_s]])
#'   Xs <- data_part$X_list[[ind_s]]
#'   crds <- data_part$crd_list[[ind_s]]
#'   Ws <- W_list[[ind_s]]
#'   result <- spBPS::BPS_pred(data = list(Y = Ys, X = Xs), coords = crds,
#'                             X_u = X_new, crd_u = crd_new,
#'                             priors = list(mu_b = matrix(rep(0, p)),
#'                                           V_b = diag(10, p),
#'                                           a = 2,
#'                                           b = 2),
#'                                           hyperpar = list(delta = delta_seq,
#'                                                           phi = phi_seq),
#'                                           W = Ws, R = 1)
#'
#'   predictions[[r]] <- result}
#'
#' }
#'
#' @export
BPS_pred <- function(data, X_u, priors, coords, crd_u, hyperpar, W, R) {
    .Call(`_spBPS_BPS_pred`, data, X_u, priors, coords, crd_u, hyperpar, W, R)
}

#' Perform the BPS sampling from posterior and posterior predictive given a set of stacking weights
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param X_u [matrix] unobserved instances covariate matrix
#' @param priors [list] priors: named \eqn{\mu_b},\eqn{V_b},\eqn{a},\eqn{b}
#' @param coords [matrix] sample coordinates for X and Y
#' @param crd_u [matrix] unboserved instances coordinates
#' @param hyperpar [list] two elemets: first named \eqn{\delta}, second named \eqn{\phi}
#' @param W [matrix] set of stacking weights
#' @param R [integer] number of desired samples
#'
#' @return [list] BPS posterior predictive samples
#'
#' @examples
#' \donttest{
#' ## Generate subsets of data
#' n <- 100
#' p <- 3
#' X <- matrix(rnorm(n*p), nrow = n, ncol = p)
#' Y <- matrix(rnorm(n), nrow = n, ncol = 1)
#' crd <- matrix(runif(n*2), nrow = n, ncol = 2)
#' data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10)
#'
#' ## Select competetive set of values for hyperparameters
#' delta_seq <- c(0.1, 0.2, 0.3)
#' phi_seq <- c(3, 4, 5)
#' ## Fit local models
#' fit_list <- vector(length = 10, mode = "list")
#' for (i in 1:10) {
#'   Yi <- data_part$Y_list[[i]]
#'   Xi <- data_part$X_list[[i]]
#'   crd_i <- data_part$crd_list[[i]]
#'   p <- ncol(Xi)
#'   bps <- spBPS::BPS_weights(data = list(Y = Yi, X = Xi),
#'                            priors = list(mu_b = matrix(rep(0, p)),
#'                                          V_b = diag(10, p),
#'                                          a = 2,
#'                                          b = 2), coords = crd_i,
#'                                          hyperpar = list(delta = delta_seq,
#'                                                          phi = phi_seq),
#'                                                          K = 5)
#'   w_hat <- bps$W
#'   epd <- bps$epd
#'   fit_list[[i]] <- list(epd, w_hat) }
#'
#' ## Model combination weights between partitions using Bayesian Predictive Stacking
#' comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1)
#' Wbps <- comb_bps$W
#' W_list <- comb_bps$W_list
#'
#' ## Generate prediction points
#' m <- 100
#' X_new <- matrix(rnorm(m*p), nrow = m, ncol = p)
#' crd_new <- matrix(runif(m*2), nrow = m, ncol = 2)
#'
#' ## Perform posterior and posterior predictive sampling
#' R <- 250
#' subset_ind <- sample(1:10, R, TRUE, Wbps)
#' postsmp_and_pred <- vector(length = R, mode = "list")
#' for (r in 1:R) {
#'  ind_s <- subset_ind[r]
#'  Ys <- matrix(data_part$Y_list[[ind_s]])
#'  Xs <- data_part$X_list[[ind_s]]
#'  crds <- data_part$crd_list[[ind_s]]
#'  Ws <- W_list[[ind_s]]
#'  result <- spBPS::BPS_post(data = list(Y = Ys, X = Xs), coords = crds,
#'                            X_u = X_new, crd_u = crd_new,
#'                            priors = list(mu_b = matrix(rep(0, p)),
#'                                          V_b = diag(10, p),
#'                                          a = 2,
#'                                          b = 2),
#'                                          hyperpar = list(delta = delta_seq,
#'                                                          phi = phi_seq),
#'                                                          W = Ws, R = 1)
#'  postsmp_and_pred[[r]] <- result}
#'
#' }
#'
#' @export
BPS_post <- function(data, X_u, priors, coords, crd_u, hyperpar, W, R) {
    .Call(`_spBPS_BPS_post`, data, X_u, priors, coords, crd_u, hyperpar, W, R)
}

#' Compute the BPS posterior samples given a set of stacking weights
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param priors [list] priors: named \eqn{\mu_b},\eqn{V_b},\eqn{a},\eqn{b}
#' @param coords [matrix] sample coordinates for X and Y
#' @param hyperpar [list] two elemets: first named \eqn{\delta}, second named \eqn{\phi}
#' @param W [matrix] set of stacking weights
#' @param R [integer] number of desired samples
#'
#' @return [matrix] BPS posterior samples
#'
BPS_postdraws <- function(data, priors, coords, hyperpar, W, R) {
    .Call(`_spBPS_BPS_postdraws`, data, priors, coords, hyperpar, W, R)
}

#' Perform prediction for BPS accelerated models - loop over prediction set
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param X_u [matrix] unobserved instances covariate matrix
#' @param priors [list] priors: named \eqn{\mu_b},\eqn{V_b},\eqn{a},\eqn{b}
#' @param coords [matrix] sample coordinates for X and Y
#' @param crd_u [matrix] unboserved instances coordinates
#' @param hyperpar [list] two elemets: first named \eqn{\delta}, second named \eqn{\phi}
#' @param W [matrix] set of stacking weights
#' @param R [integer] number of desired samples
#' @param J [integer] number of desired partition of prediction set
#'
#' @return [list] BPS posterior predictive samples
#'
spPredict_BPS <- function(data, X_u, priors, coords, crd_u, hyperpar, W, R, J) {
    .Call(`_spBPS_spPredict_BPS`, data, X_u, priors, coords, crd_u, hyperpar, W, R, J)
}

#' Compute the parameters for the posteriors distribution of \eqn{\beta} and \eqn{\Sigma} (i.e. updated parameters)
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param priors [list] priors: named \eqn{\mu_B},\eqn{V_r},\eqn{\Psi},\eqn{\nu}
#' @param coords [matrix] sample coordinates for X and Y
#' @param hyperpar [list] two elemets: first named \eqn{\alpha}, second named \eqn{\phi}
#'
#' @return [list] posterior update parameters
#'
fit_cpp_MvT <- function(data, priors, coords, hyperpar) {
    .Call(`_spBPS_fit_cpp_MvT`, data, priors, coords, hyperpar)
}

#' Sample R draws from the posterior distributions
#'
#' @param poster [list] output from \code{fit_cpp} function
#' @param R [integer] number of posterior samples
#' @param par if \code{TRUE} only \eqn{\beta} and \eqn{\Sigma} are sampled (\eqn{\omega} is omitted)
#' @param p [integer] if \code{par = TRUE}, it specifies the column number of \eqn{X}
#'
#' @return [list] posterior samples
#'
post_draws_MvT <- function(poster, R, par, p) {
    .Call(`_spBPS_post_draws_MvT`, poster, R, par, p)
}

#' Draw from the joint posterior predictive for a set of unobserved covariates
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param X_u [matrix] unobserved instances covariate matrix
#' @param d_u [matrix] unobserved instances distance matrix
#' @param d_us [matrix] cross-distance between unobserved and observed instances matrix
#' @param hyperpar [list] two elemets: first named \eqn{\alpha}, second named \eqn{\phi}
#' @param poster [list] output from \code{fit_cpp} function
#' @param R [integer] number of posterior predictive samples
#'
#' @return [list] posterior predictive samples
#'
r_pred_joint_MvT <- function(data, X_u, d_u, d_us, hyperpar, poster, R) {
    .Call(`_spBPS_r_pred_joint_MvT`, data, X_u, d_u, d_us, hyperpar, poster, R)
}

#' Draw from the joint posterior predictive for a set of unobserved covariates
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param X_u [matrix] unobserved instances covariate matrix
#' @param d_u [matrix] unobserved instances distance matrix
#' @param d_us [matrix] cross-distance between unobserved and observed instances matrix
#' @param hyperpar [list] two elemets: first named \eqn{\alpha}, second named \eqn{\phi}
#' @param poster [list] output from \code{fit_cpp} function
#' @param R [integer] number of posterior predictive samples
#'
#' @return [list] posterior predictive samples
#'
r_pred_marg_MvT <- function(data, X_u, d_u, d_us, hyperpar, poster, R) {
    .Call(`_spBPS_r_pred_marg_MvT`, data, X_u, d_u, d_us, hyperpar, poster, R)
}

#' Draw from the conditional posterior predictive for a set of unobserved covariates
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param X_u [matrix] unobserved instances covariate matrix
#' @param d_u [matrix] unobserved instances distance matrix
#' @param d_us [matrix] cross-distance between unobserved and observed instances matrix
#' @param hyperpar [list] two elemets: first named \eqn{\alpha}, second named \eqn{\phi}
#' @param poster [list] output from \code{fit_cpp_MvT} function
#' @param post [list] output from \code{post_draws_MvT} function
#'
#' @return [list] posterior predictive samples
#'
r_pred_cond_MvT <- function(data, X_u, d_u, d_us, hyperpar, poster, post) {
    .Call(`_spBPS_r_pred_cond_MvT`, data, X_u, d_u, d_us, hyperpar, poster, post)
}

#' Evaluate the density of a set of unobserved response with respect to the conditional posterior predictive
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param X_u [matrix] unobserved instances covariate matrix
#' @param Y_u [matrix] unobserved instances response matrix
#' @param d_u [matrix] unobserved instances distance matrix
#' @param d_us [matrix] cross-distance between unobserved and observed instances matrix
#' @param hyperpar [list] two elemets: first named \eqn{\alpha}, second named \eqn{\phi}
#' @param poster [list] output from \code{fit_cpp} function
#'
#' @return [double] posterior predictive density evaluation
#'
d_pred_cpp_MvT <- function(data, X_u, Y_u, d_u, d_us, hyperpar, poster) {
    .Call(`_spBPS_d_pred_cpp_MvT`, data, X_u, Y_u, d_u, d_us, hyperpar, poster)
}

#' Compute the LOOCV of the density evaluations for fixed values of the hyperparameters
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param priors [list] priors: named \eqn{\mu_B},\eqn{V_r},\eqn{\Psi},\eqn{\nu}
#' @param coords [matrix] sample coordinates for X and Y
#' @param hyperpar [list] two elemets: first named \eqn{\alpha}, second named \eqn{\phi}
#'
#' @return [vector] posterior predictive density evaluations
#'
dens_loocv_MvT <- function(data, priors, coords, hyperpar) {
    .Call(`_spBPS_dens_loocv_MvT`, data, priors, coords, hyperpar)
}

#' Compute the KCV of the density evaluations for fixed values of the hyperparameters
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param priors [list] priors: named \eqn{\mu_B},\eqn{V_r},\eqn{\Psi},\eqn{\nu}
#' @param coords [matrix] sample coordinates for X and Y
#' @param hyperpar [list] two elemets: first named \eqn{\alpha}, second named \eqn{\phi}
#' @param K [integer] number of folds
#'
#' @return [vector] posterior predictive density evaluations
#'
dens_kcv_MvT <- function(data, priors, coords, hyperpar, K) {
    .Call(`_spBPS_dens_kcv_MvT`, data, priors, coords, hyperpar, K)
}

#' Return the CV predictive density evaluations for all the model combinations
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param priors [list] priors: named \eqn{\mu_B},\eqn{V_r},\eqn{\Psi},\eqn{\nu}
#' @param coords [matrix] sample coordinates for X and Y
#' @param hyperpar [list] two elemets: first named \eqn{\alpha}, second named \eqn{\phi}
#' @param useKCV if \code{TRUE} K-fold cross validation is used instead of LOOCV (no \code{default})
#' @param K [integer] number of folds
#'
#' @return [matrix] posterior predictive density evaluations (each columns represent a different model)
#'
models_dens_MvT <- function(data, priors, coords, hyperpar, useKCV, K) {
    .Call(`_spBPS_models_dens_MvT`, data, priors, coords, hyperpar, useKCV, K)
}

#' Compute the BPS weights by convex optimization
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param priors [list] priors: named \eqn{\mu_B},\eqn{V_r},\eqn{\Psi},\eqn{\nu}
#' @param coords [matrix] sample coordinates for X and Y
#' @param hyperpar [list] two elemets: first named \eqn{\alpha}, second named \eqn{\phi}
#' @param K [integer] number of folds
#'
#' @return [matrix] posterior predictive density evaluations (each columns represent a different model)
#'
#' @examples
#' \donttest{
#' ## Generate subsets of data
#' n <- 100
#' p <- 3
#' q <- 2
#' X <- matrix(rnorm(n*p), nrow = n, ncol = p)
#' Y <- matrix(rnorm(n*q), nrow = n, ncol = q)
#' crd <- matrix(runif(n*2), nrow = n, ncol = 2)
#'
#' ## Select competitive set of values for hyperparameters
#' alfa_seq <- c(0.7, 0.8, 0.9)
#' phi_seq <- c(3, 4, 5)
#'
#' ## Perform Bayesian Predictive Stacking within subsets
#' bps <- spBPS::BPS_weights_MvT(data = list(Y = Y, X = X),
#'                               priors = list(mu_B = matrix(0, nrow = p, ncol = q),
#'                                             V_r = diag(10, p),
#'                                             Psi = diag(1, q),
#'                                             nu = 3), coords = crd,
#'                                             hyperpar = list(alpha = alfa_seq,
#'                                                             phi = phi_seq),
#'                                             K = 5)
#' }
#'
#' @export
BPS_weights_MvT <- function(data, priors, coords, hyperpar, K) {
    .Call(`_spBPS_BPS_weights_MvT`, data, priors, coords, hyperpar, K)
}

#' Compute the BPS spatial prediction given a set of stacking weights
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param X_u [matrix] unobserved instances covariate matrix
#' @param priors [list] priors: named \eqn{\mu_B},\eqn{V_r},\eqn{\Psi},\eqn{\nu}
#' @param coords [matrix] sample coordinates for X and Y
#' @param crd_u [matrix] unboserved instances coordinates
#' @param hyperpar [list] two elemets: first named \eqn{\alpha}, second named \eqn{\phi}
#' @param W [matrix] set of stacking weights
#' @param R [integer] number of desired samples
#'
#' @return [list] BPS posterior predictive samples
#'
#' @examples
#' \donttest{
#' ## Generate subsets of data
#' n <- 100
#' p <- 3
#' q <- 2
#' X <- matrix(rnorm(n*p), nrow = n, ncol = p)
#' Y <- matrix(rnorm(n*q), nrow = n, ncol = q)
#' crd <- matrix(runif(n*2), nrow = n, ncol = 2)
#' data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10)
#'
#' ## Select competitive set of values for hyperparameters
#' alfa_seq <- c(0.7, 0.8, 0.9)
#' phi_seq <- c(3, 4, 5)
#'
#' ## Fit local models
#' fit_list <- vector(length = 10, mode = "list")
#' for (i in 1:10) {
#'     Yi <- data_part$Y_list[[i]]
#'     Xi <- data_part$X_list[[i]]
#'     crd_i <- data_part$crd_list[[i]]
#'     bps <- spBPS::BPS_weights_MvT(data = list(Y = Yi, X = Xi),
#'                               priors = list(mu_B = matrix(0, nrow = p, ncol = q),
#'                                             V_r = diag(10, p),
#'                                             Psi = diag(1, q),
#'                                             nu = 3), coords = crd_i,
#'                                             hyperpar = list(alpha = alfa_seq,
#'                                                             phi = phi_seq),
#'                                             K = 5)
#'      w_hat <- bps$W
#'      epd <- bps$epd
#'      fit_list[[i]] <- list(epd, w_hat) }
#'
#' ## Model combination weights between partitions using Bayesian Predictive Stacking
#' comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1)
#' Wbps <- comb_bps$W
#' W_list <- comb_bps$W_list
#'
#' ## Generate prediction points
#' m <- 100
#' X_new <- matrix(rnorm(m*p), nrow = m, ncol = p)
#' crd_new <- matrix(runif(m*2), nrow = m, ncol = 2)
#'
#' ## Perform posterior predictive sampling
#' R <- 250
#' subset_ind <- sample(1:10, R, TRUE, Wbps)
#' predictions <- vector(length = R, mode = "list")
#' for (r in 1:R) {
#'   ind_s <- subset_ind[r]
#'   Ys <- data_part$Y_list[[ind_s]]
#'   Xs <- data_part$X_list[[ind_s]]
#'   crds <- data_part$crd_list[[ind_s]]
#'   Ws <- W_list[[ind_s]]
#'   result <- spBPS::BPS_pred_MvT(data = list(Y = Ys, X = Xs), coords = crds,
#'                                 X_u = X_new, crd_u = crd_new,
#'                                 priors = list(mu_B = matrix(0, nrow = p, ncol = q),
#'                                               V_r = diag(10, p),
#'                                               Psi = diag(1, q),
#'                                               nu = 3),
#'                                               hyperpar = list(alpha = alfa_seq,
#'                                                               phi = phi_seq),
#'                                               W = Ws, R = 1)
#'
#'   predictions[[r]] <- result}
#'
#' }
#'
#' @export
BPS_pred_MvT <- function(data, X_u, priors, coords, crd_u, hyperpar, W, R) {
    .Call(`_spBPS_BPS_pred_MvT`, data, X_u, priors, coords, crd_u, hyperpar, W, R)
}

#' Perform the BPS sampling from posterior and posterior predictive given a set of stacking weights
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param X_u [matrix] unobserved instances covariate matrix
#' @param priors [list] priors: named \eqn{\mu_B},\eqn{V_r},\eqn{\Psi},\eqn{\nu}
#' @param coords [matrix] sample coordinates for X and Y
#' @param crd_u [matrix] unboserved instances coordinates
#' @param hyperpar [list] two elemets: first named \eqn{\alpha}, second named \eqn{\phi}
#' @param W [matrix] set of stacking weights
#' @param R [integer] number of desired samples
#'
#' @return [list] BPS posterior predictive samples
#'
#' @examples
#' \donttest{
#' ## Generate subsets of data
#' n <- 100
#' p <- 3
#' q <- 2
#' X <- matrix(rnorm(n*p), nrow = n, ncol = p)
#' Y <- matrix(rnorm(n*q), nrow = n, ncol = q)
#' crd <- matrix(runif(n*2), nrow = n, ncol = 2)
#' data_part <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10)
#'
#' ## Select competitive set of values for hyperparameters
#' alfa_seq <- c(0.7, 0.8, 0.9)
#' phi_seq <- c(3, 4, 5)
#'
#' ## Fit local models
#' fit_list <- vector(length = 10, mode = "list")
#' for (i in 1:10) {
#'     Yi <- data_part$Y_list[[i]]
#'     Xi <- data_part$X_list[[i]]
#'     crd_i <- data_part$crd_list[[i]]
#'     bps <- spBPS::BPS_weights_MvT(data = list(Y = Yi, X = Xi),
#'                               priors = list(mu_B = matrix(0, nrow = p, ncol = q),
#'                                             V_r = diag(10, p),
#'                                             Psi = diag(1, q),
#'                                             nu = 3), coords = crd_i,
#'                                             hyperpar = list(alpha = alfa_seq,
#'                                                             phi = phi_seq),
#'                                             K = 5)
#'      w_hat <- bps$W
#'      epd <- bps$epd
#'      fit_list[[i]] <- list(epd, w_hat) }
#'
#' ## Model combination weights between partitions using Bayesian Predictive Stacking
#' comb_bps <- BPS_combine(fit_list = fit_list, K = 10, rp = 1)
#' Wbps <- comb_bps$W
#' W_list <- comb_bps$W_list
#'
#' ## Generate prediction points
#' m <- 100
#' X_new <- matrix(rnorm(m*p), nrow = m, ncol = p)
#' crd_new <- matrix(runif(m*2), nrow = m, ncol = 2)
#'
#' ## Perform posterior and posterior predictive sampling
#' R <- 250
#' subset_ind <- sample(1:10, R, TRUE, Wbps)
#' postsmp_and_pred <- vector(length = R, mode = "list")
#' for (r in 1:R) {
#'   ind_s <- subset_ind[r]
#'   Ys <- data_part$Y_list[[ind_s]]
#'   Xs <- data_part$X_list[[ind_s]]
#'   crds <- data_part$crd_list[[ind_s]]
#'   Ws <- W_list[[ind_s]]
#'   result <- spBPS::BPS_post_MvT(data = list(Y = Ys, X = Xs), coords = crds,
#'                                 X_u = X_new, crd_u = crd_new,
#'                                 priors = list(mu_B = matrix(0, nrow = p, ncol = q),
#'                                               V_r = diag(10, p),
#'                                               Psi = diag(1, q),
#'                                               nu = 3),
#'                                               hyperpar = list(alpha = alfa_seq,
#'                                                               phi = phi_seq),
#'                                               W = Ws, R = 1)
#'
#'   postsmp_and_pred[[r]] <- result}
#'
#' }
#'
#' @export
BPS_post_MvT <- function(data, X_u, priors, coords, crd_u, hyperpar, W, R) {
    .Call(`_spBPS_BPS_post_MvT`, data, X_u, priors, coords, crd_u, hyperpar, W, R)
}

#' Compute the BPS posterior samples given a set of stacking weights
#'
#' @param data [list] two elements: first named \eqn{Y}, second named \eqn{X}
#' @param priors [list] priors: named \eqn{\mu_B},\eqn{V_r},\eqn{\Psi},\eqn{\nu}
#' @param coords [matrix] sample coordinates for X and Y
#' @param hyperpar [list] two elemets: first named \eqn{\alpha}, second named \eqn{\phi}
#' @param W [matrix] set of stacking weights
#' @param R [integer] number of desired samples
#' @param par if \code{TRUE} only \eqn{\beta} and \eqn{\Sigma} are sampled (\eqn{\omega} is omitted)
#'
#' @return [matrix] BPS posterior samples
#'
BPS_postdraws_MvT <- function(data, priors, coords, hyperpar, W, R, par) {
    .Call(`_spBPS_BPS_postdraws_MvT`, data, priors, coords, hyperpar, W, R, par)
}

#' Compute the Euclidean distance matrix
#'
#' @param X [matrix] (tipically of \eqn{N} coordindates on \eqn{\mathbb{R}^2} )
#'
#' @return [matrix] distance matrix of the elements of \eqn{X}
#'
#' @examples
#' ## Compute the Distance matrix of dimension (n x n)
#' n <- 100
#' p <- 2
#' X <- matrix(runif(n*p), nrow = n, ncol = p)
#' distance.matrix <- arma_dist(X)
#'
#' @export
arma_dist <- function(X) {
    .Call(`_spBPS_arma_dist`, X)
}

#' Build a grid from two vector (i.e. equivalent to \code{expand.grid()} in \code{R})
#'
#' @param x [vector] first vector of numeric elements
#' @param y [vector] second vector of numeric elements
#'
#' @return [matrix] expanded grid of combinations
#'
#' @examples
#' ## Create a matrix from all combination of vectors
#' x <- seq(0, 10, length.out = 100)
#' y <- seq(-1, 1, length.out = 20)
#' grid <- expand_grid_cpp(x = x, y = y)
#'
#' @export
expand_grid_cpp <- function(x, y) {
    .Call(`_spBPS_expand_grid_cpp`, x, y)
}

#' Function to sample integers (index)
#'
#' @param size [integer] dimension of the set to sample
#' @param length [integer] number of elements to sample
#' @param p [vector] sampling probabilities
#'
#' @return [vector] sample of integers
#'
sample_index <- function(size, length, p) {
    .Call(`_spBPS_sample_index`, size, length, p)
}

#' Function to subset data for meta-analysis
#'
#' @param data [list] three elements: first named \eqn{Y}, second named \eqn{X}, third named \eqn{crd}
#' @param K [integer] number of desired subsets
#'
#' @return [list] subsets of data, and the set of indexes
#'
#' @examples
#' ## Create a list of K random subsets given a list with Y, X, and crd
#' n <- 100
#' p <- 3
#' q <- 2
#' X <- matrix(rnorm(n*p), nrow = n, ncol = p)
#' Y <- matrix(rnorm(n*q), nrow = n, ncol = q)
#' crd <- matrix(runif(n*2), nrow = n, ncol = 2)
#' subsets <- subset_data(data = list(Y = Y, X = X, crd = crd), K = 10)
#'
#' @export
subset_data <- function(data, K) {
    .Call(`_spBPS_subset_data`, data, K)
}

#' Function to subset data for meta-analysis
#'
#' @param mat [matrix] not-symmetric matrix
#'
#' @return [matrix] symmetric matrix (lower triangular of \code{mat} is used)
#'
#' @examples
#' ## Force matrix to be symmetric (avoiding numerical problems)
#' n <- 4
#' X <- matrix(runif(n*n), nrow = n, ncol = n)
#' X <- forceSymmetry_cpp(mat = X)
#'
#' @export
forceSymmetry_cpp <- function(mat) {
    .Call(`_spBPS_forceSymmetry_cpp`, mat)
}