R/RcppExports.R

In rchan26/hierarchicalFusion: Divide-and-Conquer Monte Carlo Fusion

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

construct_V_cpp <- function(s, t, end_time, C, d, precondition_matrices, Lambda) {
    .Call(`_DCFusion_construct_V_cpp`, s, t, end_time, C, d, precondition_matrices, Lambda)
}

construct_M_cpp <- function(s, t, end_time, C, d, sub_posterior_samples, sub_posterior_mean) {
    .Call(`_DCFusion_construct_M_cpp`, s, t, end_time, C, d, sub_posterior_samples, sub_posterior_mean)
}

#' Calculate the weighted mean (univariate)
#' 
#' Calculation of weighted mean when the target is univariate
#'
#' @param x vector of values
#' @param weights vector of weights
#'
#' @return the weighted mean of x
#' 
#' @examples
#' x <- c(0.4, 0.2, 0.5, 0.9, 1.4)
#' w <- c(3, 2, 5, 1, 2)
#' weighted_mean_univariate(x = x,
#'                          weights = w)
#' # returns the same using weighted.mean(x, w) function in R
#' weighted.mean(x, w)
#'
#' weighted_mean_univariate(x = x,
#'                          weights = rep(1, 5))
#' #returns the same using standard mean(x) function in R
#' mean(x)
weighted_mean_univariate <- function(x, weights) {
    .Call(`_DCFusion_weighted_mean_univariate`, x, weights)
}

#' Calculate the logarithm of rho (univariate)
#' 
#' Calculation of the log of rho acceptance probability or weight when target is
#' univariate
#'
#' @param x vector of sampled sub-posterior values
#' @param x_mean weighted mean of sampled sub-posterior values
#' @param time time T for fusion algorithm
#' @param precondition_values precondition values associated to each sub-posterior
#'
#' @return the logarithm of rho
#' 
#' @examples
#' x <- rnorm(4, 0, 1)
#' precondition_vals <- c(1, 2, 3, 4)
#' x_mean <- weighted_mean_univariate(x = x,
#'                                    weights = 1/precondition_vals)
#' log_rho_univariate(x = x,
#'                    x_mean = x_mean,
#'                    time = 0.5,
#'                    precondition_values = precondition_vals)
log_rho_univariate <- function(x, x_mean, time, precondition_values) {
    .Call(`_DCFusion_log_rho_univariate`, x, x_mean, time, precondition_values)
}

#' Calculate the variance of numbers (univariate)
#' 
#' Calculation of the weighted variance of numbers
#'
#' @param x vector of sampled sub-posterior values
#' @param x_mean weighted mean of sampled sub-posterior values
#' @param precondition_values precondition values associated to each sub-posterior
#'
#' @return the weighted variance of numbers
#'
#' @examples
#' x <- rnorm(4, 0, 1)
#' precondition_vals <- c(1, 2, 3, 4)
#' x_mean <- weighted_mean_univariate(x = x,
#'                                    weights = 1/precondition_vals)
#' weighted_variance_univariate(x = x,
#'                              x_mean = x_mean,
#'                              precondition_values = precondition_vals)
weighted_variance_univariate <- function(x, x_mean, precondition_values) {
    .Call(`_DCFusion_weighted_variance_univariate`, x, x_mean, precondition_values)
}

#' Calculate approximation to expectation of nu_j (univariate)
#' 
#' Calculation of the scaled/weighted average variation of the C trajectories
#' with respect to their individual sub-posterior means
#'
#' @param list where x_samples[[i]] ith collection of the C trajectories
#' @param sub_posterior_means vector of length C of sub-posterior means
#' @param precondition_values precondition values associated to each sub-posterior
#'
#' @return the approximated expectation of nu_j
#'
#' @examples
#' # x_samples has 5 samples and C=4
#' N <- 10
#' C <- 4
#' x_samples <- lapply(1:N, function(i) rnorm(C))
#' normalised_weights <- rep(1/N, N)
#' sub_posterior_means <- rnorm(C)
#' precond <- 1:C
#' weighted_trajectory_variation_univariate(x_samples = x_samples,
#'                                          normalised_weights = normalised_weights,
#'                                          sub_posterior_means = sub_posterior_means,
#'                                          precondition_values = precond)
#' # should be equal to the result of this:
#' sum(sapply(1:N, function(i) {
#'   sum((((x_samples[[i]]-sub_posterior_means)^2)/precond))/C
#' }))/N
weighted_trajectory_variation_univariate <- function(x_samples, normalised_weights, sub_posterior_means, precondition_values) {
    .Call(`_DCFusion_weighted_trajectory_variation_univariate`, x_samples, normalised_weights, sub_posterior_means, precondition_values)
}

compute_max_E_nu_j_univariate <- function(N, sub_posterior_samples, log_weights, time, sub_posterior_means, precondition_values) {
    .Call(`_DCFusion_compute_max_E_nu_j_univariate`, N, sub_posterior_samples, log_weights, time, sub_posterior_means, precondition_values)
}

#' Calculate the inverse of a sum of matrices
#'
#' Calculation of the inverse of a sum of a list of matrices
#'
#' @param matrices list of matrices (of same dimension)
#'
#' @return the inverse of the sum of the matrices
#'
#' @examples
#' m1 <- matrix(c(1,2,3,4), nrow = 2, ncol = 2)
#' m2 <- matrix(c(5,6,7,8), nrow = 2, ncol = 2)
#' m3 <- matrix(c(9,10,11,12), nrow = 2, ncol = 2)
#' inverse_sum_matrices(list(m1, m2, m3))
#' # returns the same as using solve() in R
#' solve(m1+m2+m3)
inverse_sum_matrices <- function(matrices) {
    .Call(`_DCFusion_inverse_sum_matrices`, matrices)
}

#' Calculate the weighted mean (multivariate)
#'
#' Calculation of weighted mean when the target is multivariate
#' 
#' @param matrix an (m x n) matrix where the ith row is the ith sample
#' @param weights list of m matrices of the same dimension (n x n)
#' @param inverse_sum_weights the inverse of the sum of the weights (can be 
#'                        calculated by passing in weights to inverse_sum_matrices)
#'
#' @return proposal mean vector
#' 
#' @examples
#' m1 <- matrix(c(1,2,3,4), nrow = 2, ncol = 2)
#' m2 <- matrix(c(5,6,7,8), nrow = 2, ncol = 2)
#' m3 <- matrix(c(9,10,11,12), nrow = 2, ncol = 2)
#' X <- matrix(c(1,2,3,4,5,6), nrow = 3, ncol = 2)
#' weighted_mean_multivariate(matrix = X,
#'                            weights = list(m1, m2, m3),
#'                            inverse_sum_weights = inverse_sum_matrices(list(m1, m2, m3)))
weighted_mean_multivariate <- function(matrix, weights, inverse_sum_weights) {
    .Call(`_DCFusion_weighted_mean_multivariate`, matrix, weights, inverse_sum_weights)
}

#' Calculate the proposal covariance matrix
#' 
#' Calculation of the proposal covariance matrix for Monte Carlo Fusion
#'
#' @param time T for fusion algorithm
#' @param weights list of m matrices of the same dimension
#'
#' @return proposal covariance matrix
#' 
#' @examples
#' m1 <- matrix(c(1,2,3,4), nrow = 2, ncol = 2)
#' m2 <- matrix(c(5,6,7,8), nrow = 2, ncol = 2)
#' m3 <- matrix(c(9,10,11,12), nrow = 2, ncol = 2)
#' calculate_proposal_cov(time = 0.5, weights = list(m1, m2, m3))
calculate_proposal_cov <- function(time, weights) {
    .Call(`_DCFusion_calculate_proposal_cov`, time, weights)
}

#' Row-wise subtraction of a vector to rows of a matrix
#' 
#' Calculates the subtraction of a vector to each row of a matrix
#'
#' @param X matrix (n x m)
#' @param vect vector of length m
#'
#' @return an (n x m) matrix, Y, where Y[i,] = X[i,]-vect
#' 
#' @examples
#' X <- matrix(c(1,2,3,4,5,6,7,8), nrow = 4, ncol = 2, byrow = T)
#' row_wise_subtraction(X = X, vect = c(1,2))
row_wise_subtraction <- function(X, vect) {
    .Call(`_DCFusion_row_wise_subtraction`, X, vect)
}

#' Calculate the logarithm of rho (multivariate)
#' 
#' Calculation of the log of rho acceptance probability or weight when target is
#' multivariate
#'
#' @param x an (m x n) matrix where the ith row is the ith sample
#' @param x_mean a vector of length n (the weighted mean of x samples)
#' @param time time T for fusion algorithm
#' @param inv_precondition_matrices list of length m of inverse 
#'                                  preconditioning matrices
#'
#' @return the logarithm of rho
#'
#' @examples
#' # set covariance matrices
#' Sig1 <- diag(2)
#' Sig2 <- matrix(c(2, 0.5, 0.5, 2), nrow = 2, ncol = 2)
#' Sig3 <- matrix(c(4, -3.2, -3.2, 4), nrow = 2, ncol = 2)
#' # sample some x values and store in the rows
#' x <- matrix(nrow = 3, ncol = 2)
#' x[1,] <- mvrnormArma(N = 1, mu = c(0, 0), Sigma = Sig1)
#' x[2,] <- mvrnormArma(N = 1, mu = c(0, 0), Sigma = Sig2)
#' x[3,] <- mvrnormArma(N = 1, mu = c(0, 0), Sigma = Sig3)
#' # calcualte precondition matrices and their inverses
#' precondition_matrices <- list(Sig1, Sig2, Sig3)
#' inv_precondition_matrices <- lapply(precondition_matrices, solve)
#' inverse_sum_weights <- inverse_sum_matrices(precondition_matrices)
#' # calculate the weighted mean where weights are the inverse precondition matrices
#' x_mean <- weighted_mean_multivariate(matrix = x,
#'                                      weights = precondition_matrices,
#'                                      inverse_sum_weights = inverse_sum_weights)
#' # calculate logarithm of rho with time T = 0.5
#' log_rho_multivariate(x = x,
#'                      x_mean = x_mean,
#'                      time = 0.5,
#'                      inv_precondition_matrices = inv_precondition_matrices)
log_rho_multivariate <- function(x, x_mean, time, inv_precondition_matrices) {
    .Call(`_DCFusion_log_rho_multivariate`, x, x_mean, time, inv_precondition_matrices)
}

#' Calculate the variance of vectors (multivariate)
#' 
#' Calculation of the weighted variance of vectors
#'
#' @param x an (m x n) matrix where the ith row is the ith sample
#' @param x_mean a vector of length n (the weighted mean of x samples)
#' @param inv_precondition_matrices list of length m of inverse 
#'                                  preconditioning matrices
#'
#' @return the weighted variance of vectors
#'
#' @examples
#' # set covariance matrices
#' Sig1 <- diag(2)
#' Sig2 <- matrix(c(2, 0.5, 0.5, 2), nrow = 2, ncol = 2)
#' Sig3 <- matrix(c(4, -3.2, -3.2, 4), nrow = 2, ncol = 2)
#' # sample some x values and store in the rows
#' x <- matrix(nrow = 3, ncol = 2)
#' x[1,] <- mvrnormArma(N = 1, mu = c(0, 0), Sigma = Sig1)
#' x[2,] <- mvrnormArma(N = 1, mu = c(0, 0), Sigma = Sig2)
#' x[3,] <- mvrnormArma(N = 1, mu = c(0, 0), Sigma = Sig3)
#' # calculate precondition matrices and their inverses
#' precondition_matrices <- list(Sig1, Sig2, Sig3)
#' inv_precondition_matrices <- lapply(precondition_matrices, solve)
#' inverse_sum_weights <- inverse_sum_matrices(precondition_matrices)
#' # calculate the weighted mean where weights are the inverse precondition matrices
#' x_mean <- weighted_mean_multivariate(matrix = x,
#'                                      weights = precondition_matrices,
#'                                      inverse_sum_weights = inverse_sum_weights)
#' weighted_variance_multivariate(x = x,
#'                                x_mean = x_mean,
#'                                inv_precondition_matrices = inv_precondition_matrices)
weighted_variance_multivariate <- function(x, x_mean, inv_precondition_matrices) {
    .Call(`_DCFusion_weighted_variance_multivariate`, x, x_mean, inv_precondition_matrices)
}

#' Calculate approximation to expectation of nu_j (multivariate)
#' 
#' Calculation of the scaled/weighted average variation of the C trajectories
#' with respect to their individual sub-posterior means
#'
#' @param list where x_samples[[i]] ith collection of the C trajectories
#' @param sub_posterior_means matrix with C rows of sub-posterior means
#' @param inv_precondition_matrices list of length m of inverse 
#'                                  preconditioning matrices
#'
#' @return the approximated expectation of nu_j
#'
#' @examples
#' N <- 10
#' C <- 4
#' d <- 3
#' x_samples <- lapply(1:N, function(i) mvrnormArma(C, rep(0,d), diag(1,d)))
#' normalised_weights <- rep(1/N, N)
#' sub_posterior_means <- mvrnormArma(C, rep(0,d), diag(1,d))
#' precond <- lapply(1:C, function(c) diag(c, d))
#' inv_precond <- lapply(precond, solve)
#' weighted_trajectory_variation_multivariate(x_samples = x_samples,
#'                                            normalised_weights = normalised_weights,
#'                                            sub_posterior_means = sub_posterior_means,
#'                                            inv_precondition_matrices = inv_precond)
#' # should be equal to the result of this:
#' sum(sapply(1:N, function(i) {
#'   sum(sapply(1:C, function(c) {
#'     diff <- x_samples[[i]][c,]-sub_posterior_means[c,]
#'     return(t(diff) %*% inv_precond[[c]] %*% diff)
#'   }))/C
#' }))/N
weighted_trajectory_variation_multivariate <- function(x_samples, normalised_weights, sub_posterior_means, inv_precondition_matrices) {
    .Call(`_DCFusion_weighted_trajectory_variation_multivariate`, x_samples, normalised_weights, sub_posterior_means, inv_precondition_matrices)
}

compute_max_E_nu_j_multivariate <- function(N, dim, sub_posterior_samples, log_weights, time, sub_posterior_means, inv_precondition_matrices, Lambda) {
    .Call(`_DCFusion_compute_max_E_nu_j_multivariate`, N, dim, sub_posterior_samples, log_weights, time, sub_posterior_means, inv_precondition_matrices, Lambda)
}

#' Calculate the logarithm of the sum of the exponentials of the arguments
#' 
#' Calculation of the log of the sum of exponential of x, but avoids computational 
#' underflow / overflow
#' 
#' @param x vector
#'
#' @return the logarithm of the sum of the exponentials of elements of x
#'         i.e. returns log(sum(exp(x)))
#'
#' @examples
#' x <- c(1000.01, 1000.02)
#' y1 <- log(sum(exp(x)))
#' print(y1)
#' y2 <- logsumexp(x)
#' print(y2) 
logsumexp <- function(x) {
    .Call(`_DCFusion_logsumexp`, x)
}

#' Calculate the ESS
#'
#' Calculates the ESS and the normalised weights from the logarithm of weihgts
#'
#' @param log_weights the logarithm of the particle weights
#'
#' @return A list with components
#' \describe{
#'  \item{log_weights}{the logarithm of the particle weights}
#'  \item{normalised_weights}{the normalised particle weights}
#'  \item{ESS}{the effective sample size of particles (the inverse of the
#'             sum of the weights squared)}
#' }
#'
#' @examples
#' particle_ESS(log_weights = c(1000.01, 1000.02, 1000.03))
particle_ESS <- function(log_weights) {
    .Call(`_DCFusion_particle_ESS`, log_weights)
}

rho_IS_univariate_ <- function(samples_to_fuse, N, m, time, precondition_values) {
    .Call(`_DCFusion_rho_IS_univariate_`, samples_to_fuse, N, m, time, precondition_values)
}

rho_IS_multivariate_ <- function(samples_to_fuse, dim, N, m, time, inv_precondition_matrices, inverse_sum_inv_precondition_matrices) {
    .Call(`_DCFusion_rho_IS_multivariate_`, samples_to_fuse, dim, N, m, time, inv_precondition_matrices, inverse_sum_inv_precondition_matrices)
}

#' Simulate from a Multivariate Gaussian Distribution
#'
#' Produces samples from the specified multivariate Gaussian distribution
#'
#' @param N the number of samples required
#' @param mu a vector of mean values
#' @param Sigma positive-definite symmetric matrix specifying the covariance of variables
#'
#' @return samples from a multivariate Gaussian distribution
#' 
#' @examples
#' samples <- mvrnormArma(N = 10000, mu = c(0, 0), Sigma = diag(2))
mvrnormArma <- function(N, mu, Sigma) {
    .Call(`_DCFusion_mvrnormArma`, N, mu, Sigma)
}

#' Simulate from a tempered Multivariate Gaussian Distribution
#'
#' Produces samples from the specified tempered multivariate Gaussian distribution
#'
#' @param N the number of samples required
#' @param mu a vector of mean values
#' @param Sigma positive-definite symmetric matrix specifying the covariance of variables
#' @param beta inverse temperature  
#'
#' @return samples from a tempered multivariate Gaussian distribution
#'
#' @examples
#' samples <- mvrnormArma(N = 10000,
#'                        mu = c(0, 0),
#'                        Sigma = diag(2))
#' tempered_samples <- mvrnormArma_tempered(N = 10000,
#'                                          mu = c(0, 0),
#'                                          Sigma = diag(2),
#'                                          beta = 2)
mvrnormArma_tempered <- function(N, mu, Sigma, beta) {
    .Call(`_DCFusion_mvrnormArma_tempered`, N, mu, Sigma, beta)
}

#' Scaled distance between two vectors
#' 
#' Calculates the scaled distance between two vectors, i.e. calculates the norm of matrix*(x-y)
#' If matrix == identity matrix, this is just the Euclidean distance
#'
#' @param x vector
#' @param y vector
#' @param matrix matrix
#'
#' @return the scaled distance between vectors x and y with matrix 
#' 
#' @examples
#' x <- c(0.3, 0.2, 0.5, 1.2)
#' y <- c(-0.5, 0.8, 1.4, 0.9)
#' scaled_distance(x, y, diag(1, 4))
#' # should equal to the Euclidean distance:
#' sqrt(0.8^2 + 0.6^2 + 0.9^2 + 0.3^2)
scaled_distance <- function(x, y, matrix) {
    .Call(`_DCFusion_scaled_distance`, x, y, matrix)
}

#' Spectral radius of a symmetric matrix
#' 
#' Calculates the spectral radius of a symmetric matrix A (the largest absolute eigenvalue)
#' 
#' @param A matrix
#' 
#' @return The spectral radius (largest absolute eigenvalue) of A
#' 
#' @examples
#' # symmetric matrix
#' # should equal 2.5
#' spectral_radius(matrix(c(2, 0.5, 0.5, 2), nrow = 2, ncol = 2))
#' # non-symmetric matrix
#' # should equal 10
#' spectral_radius(matrix(c(9, -1, 2, -2, 8, 4, 1, 1, 8), nrow = 3, ncol = 3, byrow = T))
spectral_radius <- function(A) {
    .Call(`_DCFusion_spectral_radius`, A)
}

#' Absolute eigenvalues of a matrix
#' 
#' Calculates the absolute eigenvalues of a matrix
#' 
#' @param A matrix matrix
#' 
#' @return The absolute eigenvalues of A
#' 
#' @examples
#' # symmetric matrix
#' # should equal 2.5, 1.5
#' abs_eigenvals(matrix(c(2, 0.5, 0.5, 2), nrow = 2, ncol = 2))
#' # non-symmetric matrix
#' # should equal 10, 10, 5
#' abs_eigenvals(matrix(c(9, -1, 2, -2, 8, 4, 1, 1, 8), nrow = 3, ncol = 3, byrow = T))
abs_eigenvals <- function(A) {
    .Call(`_DCFusion_abs_eigenvals`, A)
}

#' Maximal distance to hypercube
#' 
#' Calculates the maximal distance from a vector to a hypercube
#' 
#' @param beta_hat vector
#' @param hypercube vertices matrix of hypercube vertices
#' @param transform_to_X transformation matrix to X-space (original space)
#' @param transform_to_X transformation matrix to Z-space (transformed space)
#' 
#' @return The maximal distance from beta_hat to point in hypercube
maximal_distance_hypercube_to_cv <- function(beta_hat, hypercube_vertices, transform_to_X, transform_to_Z, hypercube_centre = FALSE) {
    .Call(`_DCFusion_maximal_distance_hypercube_to_cv`, beta_hat, hypercube_vertices, transform_to_X, transform_to_Z, hypercube_centre)
}

optimise_vector_product <- function(dim, x, bessel_layers, minimise) {
    .Call(`_DCFusion_optimise_vector_product`, dim, x, bessel_layers, minimise)
}

log_BLR_gradient <- function(beta, y_labels, X, X_beta, count, prior_means, prior_variances, C) {
    .Call(`_DCFusion_log_BLR_gradient`, beta, y_labels, X, X_beta, count, prior_means, prior_variances, C)
}

log_BLR_hessian <- function(X, X_beta, count, prior_variances, C) {
    .Call(`_DCFusion_log_BLR_hessian`, X, X_beta, count, prior_variances, C)
}

ea_phi_BLR_DL_vec <- function(beta, y_labels, X, count, prior_means, prior_variances, C, precondition_mat) {
    .Call(`_DCFusion_ea_phi_BLR_DL_vec`, beta, y_labels, X, count, prior_means, prior_variances, C, precondition_mat)
}

ea_phi_BLR_DL_matrix <- function(beta, y_labels, X, count, prior_means, prior_variances, C, precondition_mat) {
    .Call(`_DCFusion_ea_phi_BLR_DL_matrix`, beta, y_labels, X, count, prior_means, prior_variances, C, precondition_mat)
}

spectral_radius_BLR <- function(beta, X, count, prior_variances, C, Lambda) {
    .Call(`_DCFusion_spectral_radius_BLR`, beta, X, count, prior_variances, C, Lambda)
}

obtain_hypercube_centre_BLR <- function(bessel_layers, transform_to_X, y_labels, X, count, prior_means, prior_variances, C) {
    .Call(`_DCFusion_obtain_hypercube_centre_BLR`, bessel_layers, transform_to_X, y_labels, X, count, prior_means, prior_variances, C)
}

spectral_radius_bound_BLR_Z <- function(dim, bessel_layers, z_hat, transformed_X, count, prior_variances, C, sqrt_Lambda) {
    .Call(`_DCFusion_spectral_radius_bound_BLR_Z`, dim, bessel_layers, z_hat, transformed_X, count, prior_variances, C, sqrt_Lambda)
}

spectral_radius_global_bound_BLR_Z <- function(dim, transformed_X, count, prior_variances, C, sqrt_Lambda) {
    .Call(`_DCFusion_spectral_radius_global_bound_BLR_Z`, dim, transformed_X, count, prior_variances, C, sqrt_Lambda)
}

ea_phi_BLR_DL_bounds <- function(beta_hat, grad_log_hat, hypercube_centre_Z, dim, transformed_X, count, prior_variances, C, transform_mats, hypercube_vertices, bessel_layers, local_bounds, hypercube_centre) {
    .Call(`_DCFusion_ea_phi_BLR_DL_bounds`, beta_hat, grad_log_hat, hypercube_centre_Z, dim, transformed_X, count, prior_variances, C, transform_mats, hypercube_vertices, bessel_layers, local_bounds, hypercube_centre)
}

gamma_NB_BLR <- function(times, h, x0, y, s, t, y_labels, X, count, prior_means, prior_variances, C, precondition_mat) {
    .Call(`_DCFusion_gamma_NB_BLR`, times, h, x0, y, s, t, y_labels, X, count, prior_means, prior_variances, C, precondition_mat)
}

log_BNBR_gradient <- function(beta, y_count, X, X_beta, count, phi_rate, prior_means, prior_variances, C) {
    .Call(`_DCFusion_log_BNBR_gradient`, beta, y_count, X, X_beta, count, phi_rate, prior_means, prior_variances, C)
}

log_BNBR_hessian <- function(y_count, X, X_beta, count, phi_rate, prior_variances, C) {
    .Call(`_DCFusion_log_BNBR_hessian`, y_count, X, X_beta, count, phi_rate, prior_variances, C)
}

ea_phi_BNBR_DL_vec <- function(beta, y_count, X, count, phi_rate, prior_means, prior_variances, C, precondition_mat) {
    .Call(`_DCFusion_ea_phi_BNBR_DL_vec`, beta, y_count, X, count, phi_rate, prior_means, prior_variances, C, precondition_mat)
}

ea_phi_BNBR_DL_matrix <- function(beta, y_count, X, count, phi_rate, prior_means, prior_variances, C, precondition_mat) {
    .Call(`_DCFusion_ea_phi_BNBR_DL_matrix`, beta, y_count, X, count, phi_rate, prior_means, prior_variances, C, precondition_mat)
}

spectral_radius_BNBR <- function(y_count, beta, X, count, phi_rate, prior_variances, C, Lambda) {
    .Call(`_DCFusion_spectral_radius_BNBR`, y_count, beta, X, count, phi_rate, prior_variances, C, Lambda)
}

log_BNBR_hessian_Z <- function(dim, y_count, transformed_X, X_beta, count, phi_rate, prior_variances, C, sqrt_Lambda) {
    .Call(`_DCFusion_log_BNBR_hessian_Z`, dim, y_count, transformed_X, X_beta, count, phi_rate, prior_variances, C, sqrt_Lambda)
}

spectral_radius_BNBR_Z <- function(dim, y_count, beta, X, transformed_X, count, phi_rate, prior_variances, C, sqrt_Lambda) {
    .Call(`_DCFusion_spectral_radius_BNBR_Z`, dim, y_count, beta, X, transformed_X, count, phi_rate, prior_variances, C, sqrt_Lambda)
}

obtain_hypercube_centre_BNBR <- function(bessel_layers, transform_to_X, y_count, X, count, prior_means, prior_variances, C, phi_rate) {
    .Call(`_DCFusion_obtain_hypercube_centre_BNBR`, bessel_layers, transform_to_X, y_count, X, count, prior_means, prior_variances, C, phi_rate)
}

obtain_G_max <- function(dim, transformed_X_vec, bessel_layers, z_hat, phi_rate, method_1 = TRUE) {
    .Call(`_DCFusion_obtain_G_max`, dim, transformed_X_vec, bessel_layers, z_hat, phi_rate, method_1)
}

spectral_radius_bound_BNBR_Z <- function(dim, bessel_layers, z_hat, y_count, transformed_X, count, phi_rate, prior_variances, C, sqrt_Lambda) {
    .Call(`_DCFusion_spectral_radius_bound_BNBR_Z`, dim, bessel_layers, z_hat, y_count, transformed_X, count, phi_rate, prior_variances, C, sqrt_Lambda)
}

spectral_radius_global_bound_BNBR_Z <- function(dim, y_count, transformed_X, count, phi_rate, prior_variances, C, sqrt_Lambda) {
    .Call(`_DCFusion_spectral_radius_global_bound_BNBR_Z`, dim, y_count, transformed_X, count, phi_rate, prior_variances, C, sqrt_Lambda)
}

ea_phi_BNBR_DL_bounds <- function(beta_hat, grad_log_hat, hypercube_centre_Z, dim, y_count, transformed_X, count, phi_rate, prior_variances, C, transform_mats, hypercube_vertices, bessel_layers, local_bounds) {
    .Call(`_DCFusion_ea_phi_BNBR_DL_bounds`, beta_hat, grad_log_hat, hypercube_centre_Z, dim, y_count, transformed_X, count, phi_rate, prior_variances, C, transform_mats, hypercube_vertices, bessel_layers, local_bounds)
}

gamma_NB_BNBR <- function(times, h, x0, y, s, t, y_count, X, count, phi_rate, prior_means, prior_variances, C, precondition_mat) {
    .Call(`_DCFusion_gamma_NB_BNBR`, times, h, x0, y, s, t, y_count, X, count, phi_rate, prior_means, prior_variances, C, precondition_mat)
}

log_BRR_gradient <- function(beta, y_resp, X, X_beta, nu, sigma, prior_means, prior_variances, C) {
    .Call(`_DCFusion_log_BRR_gradient`, beta, y_resp, X, X_beta, nu, sigma, prior_means, prior_variances, C)
}

log_BRR_hessian <- function(y_resp, X, X_beta, nu, sigma, prior_variances, C) {
    .Call(`_DCFusion_log_BRR_hessian`, y_resp, X, X_beta, nu, sigma, prior_variances, C)
}

ea_phi_BRR_DL_vec <- function(beta, y_resp, X, nu, sigma, prior_means, prior_variances, C, precondition_mat) {
    .Call(`_DCFusion_ea_phi_BRR_DL_vec`, beta, y_resp, X, nu, sigma, prior_means, prior_variances, C, precondition_mat)
}

ea_phi_BRR_DL_matrix <- function(beta, y_resp, X, nu, sigma, prior_means, prior_variances, C, precondition_mat) {
    .Call(`_DCFusion_ea_phi_BRR_DL_matrix`, beta, y_resp, X, nu, sigma, prior_means, prior_variances, C, precondition_mat)
}

spectral_radius_BRR <- function(beta, y_resp, X, nu, sigma, prior_variances, C, Lambda) {
    .Call(`_DCFusion_spectral_radius_BRR`, beta, y_resp, X, nu, sigma, prior_variances, C, Lambda)
}

obtain_hypercube_centre_BRR <- function(bessel_layers, transform_to_X, y_resp, X, nu, sigma, prior_means, prior_variances, C) {
    .Call(`_DCFusion_obtain_hypercube_centre_BRR`, bessel_layers, transform_to_X, y_resp, X, nu, sigma, prior_means, prior_variances, C)
}

spectral_radius_bound_BRR_Z <- function(dim, transformed_X, nu, sigma, prior_variances, C, sqrt_Lambda) {
    .Call(`_DCFusion_spectral_radius_bound_BRR_Z`, dim, transformed_X, nu, sigma, prior_variances, C, sqrt_Lambda)
}

ea_phi_BRR_DL_bounds <- function(beta_hat, grad_log_hat, dim, transformed_X, nu, sigma, prior_variances, C, transform_mats, hypercube_vertices) {
    .Call(`_DCFusion_ea_phi_BRR_DL_bounds`, beta_hat, grad_log_hat, dim, transformed_X, nu, sigma, prior_variances, C, transform_mats, hypercube_vertices)
}

gamma_NB_BRR <- function(times, h, x0, y, s, t, y_resp, X, nu, sigma, prior_means, prior_variances, C, precondition_mat) {
    .Call(`_DCFusion_gamma_NB_BRR`, times, h, x0, y, s, t, y_resp, X, nu, sigma, prior_means, prior_variances, C, precondition_mat)
}

ea_phi_biGaussian_DL_vec <- function(x, mean_vec, sd_vec, corr, beta, precondition_mat, transform_mat) {
    .Call(`_DCFusion_ea_phi_biGaussian_DL_vec`, x, mean_vec, sd_vec, corr, beta, precondition_mat, transform_mat)
}

ea_phi_biGaussian_DL_matrix <- function(x, mean_vec, sd_vec, corr, beta, precondition_mat, transform_mat) {
    .Call(`_DCFusion_ea_phi_biGaussian_DL_matrix`, x, mean_vec, sd_vec, corr, beta, precondition_mat, transform_mat)
}

#' Obtain bounds for phi function
#'
#' Finds the lower and upper bounds of the phi function between an interval
#'
#' @param mean_vec vector of length 2 for mean
#' @param sd_vec vector of length 2 for standard deviation
#' @param corr correlation value between component 1 and component 2
#' @param beta real value
#' @param lower vector of length 2 for the lower end of interval
#' @param upper vector of length 2 for the upper end of interval
#' @param precondition_mat precondition matrix
#' @param transform_to_Z the transformation matrix to Z-space
#' @param transform_to_X the transformation matrix to X-space
#'
#' @return A list of components
#' \describe{
#'   \item{LB}{lower bound of phi}
#'   \item{UB}{upper bound of phi}
#' }
ea_phi_biGaussian_DL_bounds <- function(mean_vec, sd_vec, corr, beta, precondition_mat, transform_to_Z, transform_to_X, lower, upper) {
    .Call(`_DCFusion_ea_phi_biGaussian_DL_bounds`, mean_vec, sd_vec, corr, beta, precondition_mat, transform_to_Z, transform_to_X, lower, upper)
}

#' Obtain the global lower bound for phi function
#'
#' Finds the global bound of the phi function between a given interval
#'
#' @param x vector of length 2
#' @param mean_vec vector of length 2 for mean
#' @param sd_vec vector of length 2 for standard deviation
#' @param corr correlation value between component 1 and component 2
#' @param beta real value
#' @param precondition_mat precondition matrix
#' @param transform_mat the transformation matrix
#'
#' @return The global lower bound of phi
ea_phi_biGaussian_DL_LB <- function(mean_vec, sd_vec, corr, beta, precondition_mat) {
    .Call(`_DCFusion_ea_phi_biGaussian_DL_LB`, mean_vec, sd_vec, corr, beta, precondition_mat)
}

gamma_NB_biGaussian <- function(times, h, x0, y, s, t, mean_vec, sd_vec, corr, beta, precondition_mat, transform_mat) {
    .Call(`_DCFusion_gamma_NB_biGaussian`, times, h, x0, y, s, t, mean_vec, sd_vec, corr, beta, precondition_mat, transform_mat)
}

#' phi-function for exp(-(((x-mean)^4)*beta)/2)
#'
#' phi-function for the Exact Algorithm for exp(-(((x-mean)^4)*beta)/2)
#'
#' @param x real value
#' @param beta beta value
#' @param mean mean value
#' @param precondition precondition value
#'
#' @return real value
#'
#' @examples
#' curve(ea_phi_exp_4_DL(x, 0, 1, 1), from = -4, to = 4)
ea_phi_exp_4_DL <- function(x, mean, beta, precondition) {
    .Call(`_DCFusion_ea_phi_exp_4_DL`, x, mean, beta, precondition)
}

#' Obtain bounds for phi function
#'
#' Finds the lower and upper bounds of the phi function between a given interval
#' 
#' @param mean mean value
#' @param beta beta value
#' @param lower lower end of interval
#' @param upper upper end of interval
#' @param precondition precondition value
#'
#' @return A list of components
#' \describe{
#'   \item{LB}{lower bound of phi}
#'   \item{UB}{upper bound of phi}
#' }
#'
#' @examples
#' mu <- 0.435
#' beta <- 0.583
#' precondition <- 1.243
#' lower <- 0
#' upper <- 1.2
#' 
#' curve(ea_phi_exp_4_DL(x, mu, beta, precondition), lower, upper, ylab = 'phi')
#' abline(v=c(lower, upper))
#' abline(h=ea_phi_exp_4_DL_bounds(mean = mu,
#'                                 beta = beta,
#'                                 precondition = precondition,
#'                                 lower = lower,
#'                                 upper = upper)
#'        col = 'red', lty = 2)
ea_phi_exp_4_DL_bounds <- function(mean, beta, precondition, lower, upper) {
    .Call(`_DCFusion_ea_phi_exp_4_DL_bounds`, mean, beta, precondition, lower, upper)
}

#' Obtain the global lower bound for phi function
#'
#' Finds the global bound of the phi function between a given interval
#' 
#' @param mean mean value
#' @param beta beta value
#' @param precondition precondition value
#'
#' @return The global lower bound of phi
#'
#' @examples
#' mu <- 0.435
#' beta <- 0.583
#' precondition <- 1.243
#' lower <- 0
#' upper <- 1.6
#' 
#' curve(ea_phi_exp_4_DL(x, mu, beta, precondition), lower, upper, ylab = 'phi')
#' abline(v=c(lower, upper))
#' abline(h=ea_phi_exp_4_DL_LB(mean = mu,
#'                             beta = beta,
#'                             precondition = precondition))
#' abline(h=ea_phi_exp_4_DL_bounds(mean = mu,
#'                                 beta = beta,
#'                                 lower = lower,
#'                                 upper = upper,
#'                                 precondition = precondition), 
#'       col = 'red', lty = 2)
ea_phi_exp_4_DL_LB <- function(mean, beta, precondition) {
    .Call(`_DCFusion_ea_phi_exp_4_DL_LB`, mean, beta, precondition)
}

gamma_NB_exp_4 <- function(times, h, x0, y, s, t, mean, beta, precondition) {
    .Call(`_DCFusion_gamma_NB_exp_4`, times, h, x0, y, s, t, mean, beta, precondition)
}

dnorm_mix_tempered_unnormalised <- function(x, w, m, s, b) {
    .Call(`_DCFusion_dnorm_mix_tempered_unnormalised`, x, w, m, s, b)
}

#' phi-function for a tempered mixture Gaussian
#'
#' phi-function for the Exact Algorithm for a tempered mixture Gaussian
#'
#' @param x real value
#' @param n_comp integer number of components of mixture Gaussian
#' @param weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#' @param beta real value
#' @param precondition precondition value
#'
#' @return real value
#'
#' @examples
#' weights <- c(0.4, 0.6)
#' means <- c(-3, 6)
#' sds <- c(1, 2)
#' beta <- 1/4
#' precondition <- 1
#' curve(dnorm_mix_tempered(x,
#'                          n_comp = 2,
#'                          weights = weights,
#'                          means = means,
#'                          sds = sds,
#'                          beta = beta),
#'       -20, 40,  n = 100000, ylim = c(-0.5, 2), ylab = 'y')
#' curve(ea_phi_mixG_DL(x,
#'                      n_comp = 2,
#'                      weights = weights,
#'                      means = means,
#'                      sds = sds,
#'                      beta = beta,
#'                      precondition = precondition),
#'       add = T, n = 100000, col = 'red')
ea_phi_mixG_DL <- function(x, n_comp, weights, means, sds, beta, precondition) {
    .Call(`_DCFusion_ea_phi_mixG_DL`, x, n_comp, weights, means, sds, beta, precondition)
}

#' Obtain the global lower bound for phi function
#'
#' Finds the global bound of the phi function between a given interval
#'
#' @param n_comp integer number of components of mixture Gaussian
#' @param weights vector: weights of mixture Gaussian
#' @param means vector: means of mixture Gassuan
#' @param sds vector: st.devs of mixture Gaussian
#' @param beta real value
#' @param normalised boolean value to determine if normalisation constant is calculated
#'
#' @return The global lower bound of phi
#'
#' @examples
#' weights <- c(0.4, 0.6)
#' means <- c(-3, 6)
#' sds <- c(1, 2)
#' beta <- 1
#' precondition <- 1
#' curve(ea_phi_mixG_DL(x,
#'                      n_comp = 2,
#'                      weights = weights,
#'                      means = means,
#'                      sds = sds,
#'                      beta = beta,
#'                      precondition = precondition),
#'       from = -20, to = 25, n = 10000,
#'       ylab = 'phi')
#' PHI <- ea_phi_mixG_DL_LB(n_comp = 2,
#'                          weights = weights,
#'                          means = means,
#'                          sds = sds,
#'                          beta = beta,
#'                          precondition = precondition,
#'                          bounds_multiplier = 1)
#' abline(h = PHI, col = 'red')
ea_phi_mixG_DL_LB <- function(n_comp, weights, means, sds, beta, precondition, bounds_multiplier) {
    .Call(`_DCFusion_ea_phi_mixG_DL_LB`, n_comp, weights, means, sds, beta, precondition, bounds_multiplier)
}

ea_phi_multiGaussian_DL_vec <- function(x, mu, inv_Sigma, beta, precondition_mat) {
    .Call(`_DCFusion_ea_phi_multiGaussian_DL_vec`, x, mu, inv_Sigma, beta, precondition_mat)
}

ea_phi_multiGaussian_DL_matrix <- function(x, mu, inv_Sigma, beta, precondition_mat) {
    .Call(`_DCFusion_ea_phi_multiGaussian_DL_matrix`, x, mu, inv_Sigma, beta, precondition_mat)
}

#' Obtain bounds for phi function
#'
#' Finds the lower and upper bounds of the phi function between an interval
#'
#' @param mu vector of length dim for mean
#' @param inv_Sigma dim x dim inverse covariance matrix
#' @param beta real value
#' @param precondition_mat dim x dim precondition matrix
#' @param hypercube_vertices matrix which determines the points to evaluate
#'                           phi which give the bounds of phi
#'
#' @return A list of components
#' \describe{
#'   \item{LB}{lower bound of phi}
#'   \item{UB}{upper bound of phi}
#' }
ea_phi_multiGaussian_DL_bounds <- function(mu, inv_Sigma, beta, precondition_mat, V) {
    .Call(`_DCFusion_ea_phi_multiGaussian_DL_bounds`, mu, inv_Sigma, beta, precondition_mat, V)
}

POE_multiGaussian_DL <- function(bessel_layers, mean, dim) {
    .Call(`_DCFusion_POE_multiGaussian_DL`, bessel_layers, mean, dim)
}

#' Obtain the global lower bound for phi function
#'
#' Finds the global bound of the phi function between a given interval
#'
#' @param mu vector of length dim for mean
#' @param inv_Sigma dim x dim inverse covariance matrix
#' @param beta real value
#' @param precondition_mat dim x dim precondition matrix
#'
#' @return The global lower bound of phi
ea_phi_multiGaussian_DL_LB <- function(mu, inv_Sigma, beta, precondition_mat) {
    .Call(`_DCFusion_ea_phi_multiGaussian_DL_LB`, mu, inv_Sigma, beta, precondition_mat)
}

gamma_NB_multiGaussian <- function(times, h, x0, y, s, t, dim, mu, inv_Sigma, beta, precondition_mat) {
    .Call(`_DCFusion_gamma_NB_multiGaussian`, times, h, x0, y, s, t, dim, mu, inv_Sigma, beta, precondition_mat)
}

#' phi-function for tempered Gaussian distribution
#'
#' phi-function for the Exact Algorithm for tempered Gaussian distribution
#'
#' @param x real value
#' @param mean mean value
#' @param sd standard deviation value
#' @param beta real value
#' @param precondition precondition value
#'
#' @return real value
#'
#' @examples
#' curve(ea_phi_uniGaussian_DL(x, 0, 1, 1, 1), from = -4, to = 4)
ea_phi_uniGaussian_DL <- function(x, mean, sd, beta, precondition) {
    .Call(`_DCFusion_ea_phi_uniGaussian_DL`, x, mean, sd, beta, precondition)
}

#' Obtain bounds for phi function
#'
#' Finds the lower and upper bounds of the phi function between an interval
#'
#' @param mean mean value
#' @param sd standard deviation value
#' @param beta real value
#' @param lower lower end of interval
#' @param upper upper end of interval
#' @param precondition precondition value
#'
#' @return A list of components
#' \describe{
#'   \item{LB}{lower bound of phi}
#'   \item{UB}{upper bound of phi}
#' }
#'
#' @examples
#' mu <- 0.423
#' sd <- 3.231
#' beta <- 0.8693
#' precondition <- 1.243
#' lower <- -2.823
#' upper <- 4.322
#' curve(ea_phi_uniGaussian_DL(x, mu, sd, beta, precondition), lower, upper)
#' abline(h=ea_phi_uniGaussian_DL_bounds(mean = mu,
#'                                       sd = sd,
#'                                       beta = beta,
#'                                       lower = lower,
#'                                       upper = upper,
#'                                       precondition = precondition),
#'        col = 'red', lty = 2)
#' 
#' # another example where the mean is not in the interval
#' mu <- 0.423
#' sd <- 3.231
#' beta <- 0.8693
#' precondition <- 1.243
#' lower <- 2.823
#' upper <- 5.322
#' curve(ea_phi_uniGaussian_DL(x, mu, sd, beta, precondition), lower, upper)
#' abline(h=ea_phi_uniGaussian_DL_bounds(mean = mu,
#'                                       sd = sd,
#'                                       beta = beta,
#'                                       precondition = precondition,
#'                                       lower = lower,
#'                                       upper = upper),
#'        col = 'red', lty = 2)
ea_phi_uniGaussian_DL_bounds <- function(mean, sd, beta, precondition, lower, upper) {
    .Call(`_DCFusion_ea_phi_uniGaussian_DL_bounds`, mean, sd, beta, precondition, lower, upper)
}

#' Obtain the global lower bound for phi function
#'
#' Finds the global bound of the phi function between a given interval
#'
#' @param mean mean value
#' @param sd standard deviation value
#' @param beta real value
#' @param precondition precondition value
#'
#' @return The global lower bound of phi
#'
#' @examples
#' mu <- 0.423
#' sd <- 3.231
#' beta <- 0.8693
#' precondition <- 1.243
#' lower <- -2.823
#' upper <- 4.322
#' curve(ea_phi_uniGaussian_DL(x, mu, sd, beta, precondition), lower, upper)
#' abline(h=ea_phi_uniGaussian_DL_bounds(mean = mu,
#'                                       sd = sd,
#'                                       beta = beta,
#'                                       precondition = precondition,
#'                                       lower = lower,
#'                                       upper = upper),
#'        col = 'red', lty = 2)
#' abline(h=ea_phi_uniGaussian_DL_LB(mean = mu,
#'                                   sd = sd,
#'                                   beta = beta, 
#'                                   precondition = precondition),
#'        col = 'blue', lty = 3)
#' 
#' # another example where the mean is not in the interval
#' mu <- 0.423
#' sd <- 3.231
#' beta <- 0.8693
#' precondition <- 1.243
#' lower <- 2.823
#' upper <- 5.322
#' curve(ea_phi_uniGaussian_DL(x, mu, sd, beta, precondition), lower, upper)
#' abline(h=ea_phi_uniGaussian_DL_bounds(mean = mu,
#'                                       sd = sd,
#'                                       beta = beta,
#'                                       lower = lower,
#'                                       upper = upper,
#'                                       precondition = precondition),
#'        col = 'red', lty = 2)
#' abline(h=ea_phi_uniGaussian_DL_LB(mean = mu,
#'                                   sd = sd,
#'                                   beta = beta, 
#'                                   precondition = precondition),
#'        col = 'blue', lty = 3)
ea_phi_uniGaussian_DL_LB <- function(mean, sd, beta, precondition) {
    .Call(`_DCFusion_ea_phi_uniGaussian_DL_LB`, mean, sd, beta, precondition)
}

gamma_NB_uniGaussian <- function(times, h, x0, y, s, t, mean, sd, beta, precondition) {
    .Call(`_DCFusion_gamma_NB_uniGaussian`, times, h, x0, y, s, t, mean, sd, beta, precondition)
}