R/RcppExports.R

In jsilve24/mongrel: Bayesian Multinomial Logistic Normal Regression

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

#' Solve Bayesian Multivariate Conjugate Linear Model
#' 
#' See details for model.  Notation: \code{N} is number of samples,
#' \code{D} is the dimension of the response, \code{Q} is number
#' of covariates. 
#' 
#' @param Y matrix of dimension D x N
#' @param X matrix of covariates of dimension Q x N
#' @param Theta matrix of prior mean of dimension D x Q
#' @param Gamma covariance matrix of dimension Q x Q
#' @param Xi covariance matrix of dimension D x D
#' @param upsilon scalar (must be > D-1) degrees of freedom for InvWishart prior
#' @param n_samples number of samples to draw (default: 2000)
#' 
#' @details 
#'    \deqn{Y ~ MN_{D-1 x N}(Lambda*X, Sigma, I_N)}
#'    \deqn{Lambda ~ MN_{D-1 x Q}(Theta, Sigma, Gamma)}
#'    \deqn{Sigma ~ InvWish(upsilon, Xi)}
#' This function provides a means of sampling from the posterior distribution of 
#' \code{Lambda} and \code{Sigma}. 
#' @return List with components 
#' 1. Lambda Array of dimension (D-1) x Q x n_samples (posterior samples)
#' 2. Sigma Array of dimension (D-1) x (D-1) x n_samples (posterior samples)
#' @export
#' @md
#' @examples
#' sim <- pibble_sim()
#' eta.hat <- t(driver::alr(t(sim$Y+0.65)))
#' fit <- conjugateLinearModel(eta.hat, sim$X, sim$Theta, sim$Gamma, 
#'                             sim$Xi, sim$upsilon, n_samples=2000)
conjugateLinearModel <- function(Y, X, Theta, Gamma, Xi, upsilon, n_samples = 2000L) {
    .Call('_stray_conjugateLinearModel', PACKAGE = 'stray', Y, X, Theta, Gamma, Xi, upsilon, n_samples)
}

#' Calculations for the Collapsed Maltipoo Model
#'
#' Functions providing access to the Log Likelihood, Gradient, and Hessian
#' of the collapsed maltipoo model. Note: These are convenience functions
#' but are not as optimized as direct coding of the MaltipooCollapsed
#' C++ class due to a lack of Memoization. By contrast function optimMaltipooCollapsed
#' is much more optimized and massively cuts down on repeated calculations.
#' A more efficient Rcpp module based implementation of these functions
#' may following if the future. For model details see \code{\link{optimMaltipooCollapsed}}
#' documentation
#' @inheritParams optimMaltipooCollapsed
#' @param eta matrix (D-1)xN of parameter values at which to calculate quantities
#' @param sylv (default:false) if true and if N < D-1 will use sylvester determinant
#'   identity to speed computation
#' @param ell P-vector of scale factors for each variance component (aka VCScale) 
#' @name loglikMaltipooCollapsed
#' @export
loglikMaltipooCollapsed <- function(Y, upsilon, Theta, X, KInv, U, eta, ell, sylv = FALSE) {
    .Call('_stray_loglikMaltipooCollapsed', PACKAGE = 'stray', Y, upsilon, Theta, X, KInv, U, eta, ell, sylv)
}

#' @rdname loglikMaltipooCollapsed
#' @export
gradMaltipooCollapsed <- function(Y, upsilon, Theta, X, KInv, U, eta, ell, sylv = FALSE) {
    .Call('_stray_gradMaltipooCollapsed', PACKAGE = 'stray', Y, upsilon, Theta, X, KInv, U, eta, ell, sylv)
}

#' @rdname loglikMaltipooCollapsed
#' @export
hessMaltipooCollapsed <- function(Y, upsilon, Theta, X, KInv, U, eta, ell, sylv = FALSE) {
    .Call('_stray_hessMaltipooCollapsed', PACKAGE = 'stray', Y, upsilon, Theta, X, KInv, U, eta, ell, sylv)
}

#' Function to Optimize the Collapsed Maltipoo Model
#' 
#' See details for model. Should likely be followed by function 
#' \code{\link{uncollapsePibble}}. Notation: \code{N} is number of samples,
#' \code{D} is number of multinomial categories, and \code{Q} is number
#' of covariates. 
#' 
#' @param Y D x N matrix of counts
#' @param upsilon (must be > D)
#' @param Theta D-1 x Q matrix the prior mean for regression coefficients
#' @param X Q x N matrix of covariates
#' @param KInv D-1 x D-1 symmetric positive-definite matrix
#' @param U a PQxQ matrix of stacked variance components 
#' @param init D-1 x N matrix of initial guess for eta used for optimization
#' @param ellinit P vector of initial guess for ell used for optimization
#' @param n_samples number of samples for Laplace Approximation (=0 very fast
#'    as no inversion or decomposition of Hessian is required)
#' @param calcGradHess if n_samples=0 should Gradient and Hessian 
#'   still be calculated using closed form solutions?
#' @param b1 (ADAM) 1st moment decay parameter (recommend 0.9) "aka momentum"
#' @param b2 (ADAM) 2nd moment decay parameter (recommend 0.99 or 0.999)
#' @param step_size (ADAM) step size for descent (recommend 0.001-0.003)
#' @param epsilon (ADAM) parameter to avoid divide by zero
#' @param eps_f (ADAM) normalized function improvement stopping criteria 
#' @param eps_g (ADAM) normalized gradient magnitude stopping criteria
#' @param max_iter (ADAM) maximum number of iterations before stopping
#' @param verbose (ADAM) if true will print stats for stopping criteria and 
#'   iteration number
#' @param verbose_rate (ADAM) rate to print verbose stats to screen
#' @param decomp_method decomposition of hessian for Laplace approximation
#'   'eigen' (more stable-slightly, slower) or 'cholesky' (less stable, faster, default)
#' @param eigvalthresh threshold for negative eigenvalues in 
#'   decomposition of negative inverse hessian (should be <=0)
#' @param jitter (default: 0) if >0 then adds that factor to diagonal of Hessian 
#' before decomposition (to improve matrix conditioning)
#'   
#' @details Notation: Let Z_j denote the J-th row of a matrix Z.
#' Model:
#'    \deqn{Y_j \sim Multinomial(Pi_j)}
#'    \deqn{Pi_j = Phi^{-1}(Eta_j)}
#'    \deqn{Eta \sim T_{D-1, N}(upsilon, Theta*X, K, A)}
#'    
#'  Where A = (I_N + e^{ell_1}*X*U_1*X' + ... + e^{ell_P}*X*U_P*X' ),
#'  K is a D-1xD-1 covariance and Phi is ALRInv_D transform. 
#' 
#' Gradient and Hessian calculations are fast as they are computed using closed
#' form solutions. That said, the Hessian matrix can be quite large 
#' \[N*(D-1) x N*(D-1)\] and storage may be an issue. 
#' 
#' Note: Warnings about large negative eigenvalues can either signal 
#' that the optimizer did not reach an optima or (more commonly in my experience)
#' that the prior / degrees of freedom for the covariance (given by parameters
#' \code{upsilon} and \code{KInv}) were too specific and at odds with the observed data.
#' If you get this warning try the following. 
#' 1. Try restarting the optimization using a different initial guess for eta
#' 2. Try decreasing (or even increasing)\code{step_size} (by increments of 0.001 or 0.002) 
#'   and increasing \code{max_iter} parameters in optimizer. Also can try 
#'   increasing \code{b1} to 0.99 and decreasing \code{eps_f} by a few orders
#'   of magnitude
#' 3. Try relaxing prior assumptions regarding covariance matrix. (e.g., may want
#' to consider decreasing parameter \code{upsilon} closer to a minimum value of 
#' D)
#' 4. Try adding small amount of jitter (e.g., set \code{jitter=1e-5}) to address
#'   potential floating point errors. 
#' @return List containing (all with respect to found optima)
#' 1. LogLik - Log Likelihood of collapsed model (up to proportionality constant)
#' 2. Gradient - (if \code{calcGradHess}=true)
#' 3. Hessian - (if \code{calcGradHess}=true) of the POSITIVE log posterior
#' 4. Pars - Parameter value of eta 
#' 5. Samples - (D-1) x N x n_samples array containing posterior samples of eta 
#'   based on Laplace approximation (if n_samples>0)
#' 6. VCScale - value of e^ell_i at optima
#' 7. logInvNegHessDet - the log determinant of the covariacne of the Laplace 
#'    approximation, useful for calculating marginal likelihood 
#' @md 
#' @export
#' @name optimMaltipooCollapsed
#' @references S. Ruder (2016) \emph{An overview of gradient descent 
#' optimization algorithms}. arXiv 1609.04747
#' @seealso \code{\link{uncollapsePibble}}
optimMaltipooCollapsed <- function(Y, upsilon, Theta, X, KInv, U, init, ellinit, n_samples = 2000L, calcGradHess = TRUE, b1 = 0.9, b2 = 0.99, step_size = 0.003, epsilon = 10e-7, eps_f = 1e-10, eps_g = 1e-4, max_iter = 10000L, verbose = FALSE, verbose_rate = 10L, decomp_method = "cholesky", eigvalthresh = 0, jitter = 0) {
    .Call('_stray_optimMaltipooCollapsed', PACKAGE = 'stray', Y, upsilon, Theta, X, KInv, U, init, ellinit, n_samples, calcGradHess, b1, b2, step_size, epsilon, eps_f, eps_g, max_iter, verbose, verbose_rate, decomp_method, eigvalthresh, jitter)
}

#' Calculations for the Collapsed Pibble Model
#'
#' Functions providing access to the Log Likelihood, Gradient, and Hessian
#' of the collapsed pibble model. Note: These are convenience functions
#' but are not as optimized as direct coding of the PibbleCollapsed
#' C++ class due to a lack of Memoization. By contrast function optimPibbleCollapsed
#' is much more optimized and massively cuts down on repeated calculations.
#' A more efficient Rcpp module based implementation of these functions
#' may following if the future. For model details see \code{\link{optimPibbleCollapsed}}
#' documentation
#'
#' @inheritParams optimPibbleCollapsed
#' @param AInv Inverse of A for LTP (for Pibble defined as 
#'   AInv = solve(diag(N)+ t(X)\%*\%Gamma\%*\%X) )
#' @param KInv Inverse of K for LTP (for Pibble defined as KInv = solve(Xi))
#' @param eta matrix (D-1)xN of parameter values at which to calculate quantities
#' @param sylv (default:false) if true and if N < D-1 will use sylvester determinant
#'   identity to speed computation
#' @return see below
#'   \itemize{
#'     \item loglikPibbleCollapsed - double
#'     \item gradPibbleCollapsed - vector
#'     \item hessPibbleCollapsed- matrix
#'   }
#' @md
#' @export
#' @examples
#' D <- 10
#' Q <- 2
#' N <- 30
#'
#' # Simulate Data
#' Sigma <- diag(sample(1:8, D-1, replace=TRUE))
#' Sigma[2, 3] <- Sigma[3,2] <- -1
#' Gamma <- diag(sqrt(rnorm(Q)^2))
#' Theta <- matrix(0, D-1, Q)
#' Phi <-  Theta + t(chol(Sigma))%*%matrix(rnorm(Q*(D-1)), nrow=D-1)%*%chol(Gamma)
#' X <- matrix(rnorm(N*(Q-1)), Q-1, N)
#' X <- rbind(1, X)
#' Eta <- Phi%*%X + t(chol(Sigma))%*%matrix(rnorm(N*(D-1)), nrow=D-1)
#' Pi <- t(driver::alrInv(t(Eta)))
#' Y <- matrix(0, D, N)
#' for (i in 1:N) Y[,i] <- rmultinom(1, sample(5000:10000), prob = Pi[,i])
#'
#' # Priors
#' upsilon <- D+10
#' Xi <- Sigma*(upsilon-D)
#'
#' # Precompute
#' KInv <- solve(Xi)
#' AInv <- solve(diag(N)+ t(X)%*%Gamma%*%X)
#' ThetaX <- Theta%*%X
#'
#'
#' loglikPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, Eta)
#' gradPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, Eta)[1:5]
#' hessPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, Eta)[1:5,1:5]
loglikPibbleCollapsed <- function(Y, upsilon, ThetaX, KInv, AInv, eta, sylv = FALSE) {
    .Call('_stray_loglikPibbleCollapsed', PACKAGE = 'stray', Y, upsilon, ThetaX, KInv, AInv, eta, sylv)
}

#' @rdname loglikPibbleCollapsed
#' @export
gradPibbleCollapsed <- function(Y, upsilon, ThetaX, KInv, AInv, eta, sylv = FALSE) {
    .Call('_stray_gradPibbleCollapsed', PACKAGE = 'stray', Y, upsilon, ThetaX, KInv, AInv, eta, sylv)
}

#' @rdname loglikPibbleCollapsed
#' @export
hessPibbleCollapsed <- function(Y, upsilon, ThetaX, KInv, AInv, eta, sylv = FALSE) {
    .Call('_stray_hessPibbleCollapsed', PACKAGE = 'stray', Y, upsilon, ThetaX, KInv, AInv, eta, sylv)
}

#' Function to Optimize the Collapsed Pibble Model
#' 
#' See details for model. Should likely be followed by function 
#' \code{\link{uncollapsePibble}}. Notation: \code{N} is number of samples,
#' \code{D} is number of multinomial categories, and \code{Q} is number
#' of covariates. 
#' 
#' @param Y D x N matrix of counts
#' @param upsilon (must be > D)
#' @param ThetaX D-1 x N matrix formed by Theta*X (Theta is Prior mean 
#'    for regression coefficients) 
#' @param KInv D-1 x D-1 precision matrix (inverse of Xi)
#' @param AInv N x N precision matrix given by (I_N + X'*Gamma*X)^{-1}
#' @param init D-1 x N matrix of initial guess for eta used for optimization
#' @param n_samples number of samples for Laplace Approximation (=0 very fast
#'    as no inversion or decomposition of Hessian is required)
#' @param calcGradHess if n_samples=0 should Gradient and Hessian 
#'   still be calculated using closed form solutions?
#' @param b1 (ADAM) 1st moment decay parameter (recommend 0.9) "aka momentum"
#' @param b2 (ADAM) 2nd moment decay parameter (recommend 0.99 or 0.999)
#' @param step_size (ADAM) step size for descent (recommend 0.001-0.003)
#' @param epsilon (ADAM) parameter to avoid divide by zero
#' @param eps_f (ADAM) normalized function improvement stopping criteria 
#' @param eps_g (ADAM) normalized gradient magnitude stopping criteria
#' @param max_iter (ADAM) maximum number of iterations before stopping
#' @param verbose (ADAM) if true will print stats for stopping criteria and 
#'   iteration number
#' @param verbose_rate (ADAM) rate to print verbose stats to screen
#' @param decomp_method decomposition of hessian for Laplace approximation
#'   'eigen' (more stable-slightly, slower) or 'cholesky' (less stable, faster, default)
#' @param optim_method (default:"adam") or "lbfgs"
#' @param eigvalthresh threshold for negative eigenvalues in 
#'   decomposition of negative inverse hessian (should be <=0)
#' @param jitter (default: 0) if >=0 then adds that factor to diagonal of Hessian 
#' before decomposition (to improve matrix conditioning)
#' @param multDirichletBoot if >0 (overrides laplace approximation) and samples
#'  eta efficiently at MAP estimate from pseudo Multinomial-Dirichlet posterior.
#' @param useSylv (default: true) if N<D-1 uses Sylvester Determinant Identity
#'   to speed up calculation of log-likelihood and gradients. 
#' @param ncores (default:-1) number of cores to use, if ncores==-1 then 
#' uses default from OpenMP typically to use all available cores. 
#' @param seed (random seed for Laplace approximation -- integer)
#'  
#' @details Notation: Let Z_j denote the J-th row of a matrix Z.
#' Model:
#'    \deqn{Y_j \sim Multinomial(Pi_j)}
#'    \deqn{Pi_j = Phi^{-1}(Eta_j)}
#'    \deqn{Eta \sim T_{D-1, N}(upsilon, Theta*X, K, A)}
#' Where A = I_N + X * Gamma * X', K is a (D-1)x(D-1) covariance 
#' matrix, Gamma is a Q x Q covariance matrix, and Phi^{-1} is ALRInv_D 
#' transform. 
#' 
#' Gradient and Hessian calculations are fast as they are computed using closed
#' form solutions. That said, the Hessian matrix can be quite large 
#' \[N*(D-1) x N*(D-1)\] and storage may be an issue. 
#' 
#' Note: Warnings about large negative eigenvalues can either signal 
#' that the optimizer did not reach an optima or (more commonly in my experience)
#' that the prior / degrees of freedom for the covariance (given by parameters
#' \code{upsilon} and \code{KInv}) were too specific and at odds with the observed data.
#' If you get this warning try the following. 
#' 1. Try restarting the optimization using a different initial guess for eta
#' 2. Try decreasing (or even increasing )\code{step_size} (by increments of 0.001 or 0.002) 
#'   and increasing \code{max_iter} parameters in optimizer. Also can try 
#'   increasing \code{b1} to 0.99 and decreasing \code{eps_f} by a few orders
#'   of magnitude
#' 3. Try relaxing prior assumptions regarding covariance matrix. (e.g., may want
#' to consider decreasing parameter \code{upsilon} closer to a minimum value of 
#' D)
#' 4. Try adding small amount of jitter (e.g., set \code{jitter=1e-5}) to address
#'   potential floating point errors. 
#' @return List containing (all with respect to found optima)
#' 1. LogLik - Log Likelihood of collapsed model (up to proportionality constant)
#' 2. Gradient - (if \code{calcGradHess}=true)
#' 3. Hessian - (if \code{calcGradHess}=true) of the POSITIVE LOG POSTERIOR
#' 4. Pars - Parameter value of eta at optima
#' 5. Samples - (D-1) x N x n_samples array containing posterior samples of eta 
#'   based on Laplace approximation (if n_samples>0)
#' 6. Timer - Vector of Execution Times
#' 7. logInvNegHessDet - the log determinant of the covariacne of the Laplace 
#'    approximation, useful for calculating marginal likelihood 
#' @md 
#' @export
#' @name optimPibbleCollapsed
#' @references S. Ruder (2016) \emph{An overview of gradient descent 
#' optimization algorithms}. arXiv 1609.04747
#' 
#' JD Silverman K Roche, ZC Holmes, LA David, S Mukherjee. 
#'   \emph{Bayesian Multinomial Logistic Normal Models through Marginally Latent Matrix-T Processes}. 
#'   2019, arXiv e-prints, arXiv:1903.11695
#' @seealso \code{\link{uncollapsePibble}}
#' @examples
#' sim <- pibble_sim()
#' 
#' # Fit model for eta
#' fit <- optimPibbleCollapsed(sim$Y, sim$upsilon, sim$Theta%*%sim$X, sim$KInv, 
#'                              sim$AInv, random_pibble_init(sim$Y))  
optimPibbleCollapsed <- function(Y, upsilon, ThetaX, KInv, AInv, init, n_samples = 2000L, calcGradHess = TRUE, b1 = 0.9, b2 = 0.99, step_size = 0.003, epsilon = 10e-7, eps_f = 1e-10, eps_g = 1e-4, max_iter = 10000L, verbose = FALSE, verbose_rate = 10L, decomp_method = "cholesky", optim_method = "adam", eigvalthresh = 0, jitter = 0, multDirichletBoot = -1.0, useSylv = TRUE, ncores = -1L, seed = -1L) {
    .Call('_stray_optimPibbleCollapsed', PACKAGE = 'stray', Y, upsilon, ThetaX, KInv, AInv, init, n_samples, calcGradHess, b1, b2, step_size, epsilon, eps_f, eps_g, max_iter, verbose, verbose_rate, decomp_method, optim_method, eigvalthresh, jitter, multDirichletBoot, useSylv, ncores, seed)
}

#' Uncollapse output from optimPibbleCollapsed to full pibble Model
#' 
#' See details for model. Should likely be called following 
#' \code{\link{optimPibbleCollapsed}}. Notation: \code{N} is number of samples,
#' \code{D} is number of multinomial categories, \code{Q} is number
#' of covariates, \code{iter} is the number of samples of \code{eta} (e.g., 
#' the parameter \code{n_samples} in the function \code{optimPibbleCollapsed})
#' 
#' @param eta array of dimension (D-1) x N x iter (e.g., \code{Pars} output of 
#'   function optimPibbleCollapsed)
#' @param X matrix of covariates of dimension Q x N
#' @param Theta matrix of prior mean of dimension (D-1) x Q
#' @param Gamma covariance matrix of dimension Q x Q
#' @param Xi covariance matrix of dimension (D-1) x (D-1)
#' @param upsilon scalar (must be > D) degrees of freedom for InvWishart prior
#' @param ret_mean if true then uses posterior mean of Lambda and Sigma 
#'   corresponding to each sample of eta rather than sampling from 
#'   posterior of Lambda and Sigma (useful if Laplace approximation
#'   is not used (or fails) in optimPibbleCollapsed)
#' @param seed seed to use for random number generation 
#' @param ncores (default:-1) number of cores to use, if ncores==-1 then 
#' uses default from OpenMP typically to use all available cores. 
#'  
#' @details Notation: Let Z_j denote the J-th row of a matrix Z.
#' While the collapsed model is given by:
#'    \deqn{Y_j ~ Multinomial(Pi_j)}
#'    \deqn{Pi_j = Phi^{-1}(Eta_j)}
#'    \deqn{Eta ~ T_{D-1, N}(upsilon, Theta*X, K, A)}
#' Where A = I_N + X * Gamma * X', K = Xi is a (D-1)x(D-1) covariance 
#' matrix, Gamma is a Q x Q covariance matrix, and Phi^{-1} is ALRInv_D 
#' transform. 
#' 
#' The uncollapsed model (Full pibble model) is given by:
#'    \deqn{Y_j ~ Multinomial(Pi_j)}
#'    \deqn{Pi_j = Phi^{-1}(Eta_j)}
#'    \deqn{Eta ~ MN_{D-1 x N}(Lambda*X, Sigma, I_N)}
#'    \deqn{Lambda ~ MN_{D-1 x Q}(Theta, Sigma, Gamma)}
#'    \deqn{Sigma ~ InvWish(upsilon, Xi)}
#' This function provides a means of sampling from the posterior distribution of 
#' \code{Lambda} and \code{Sigma} given posterior samples of \code{Eta} from 
#' the collapsed model. 
#' @return List with components 
#' 1. Lambda Array of dimension (D-1) x Q x iter (posterior samples)
#' 2. Sigma Array of dimension (D-1) x (D-1) x iter (posterior samples)
#' 3. Timer
#' @export
#' @md
#' @seealso \code{\link{optimPibbleCollapsed}}
#' @references JD Silverman K Roche, ZC Holmes, LA David, S Mukherjee. 
#'   Bayesian Multinomial Logistic Normal Models through Marginally Latent Matrix-T Processes. 
#'   2019, arXiv e-prints, arXiv:1903.11695
#' @examples
#' sim <- pibble_sim()
#' 
#' # Fit model for eta
#' fit <- optimPibbleCollapsed(sim$Y, sim$upsilon, sim$Theta%*%sim$X, sim$KInv, 
#'                              sim$AInv, random_pibble_init(sim$Y))  
#' 
#' # Finally obtain samples from Lambda and Sigma
#' fit2 <- uncollapsePibble(fit$Samples, sim$X, sim$Theta, 
#'                                    sim$Gamma, sim$Xi, sim$upsilon, 
#'                                    seed=2849)
uncollapsePibble <- function(eta, X, Theta, Gamma, Xi, upsilon, seed, ret_mean = FALSE, ncores = -1L) {
    .Call('_stray_uncollapsePibble', PACKAGE = 'stray', eta, X, Theta, Gamma, Xi, upsilon, seed, ret_mean, ncores)
}

rMatNormalCholesky_test <- function(M, LU, LV, discard) {
    .Call('_stray_rMatNormalCholesky_test', PACKAGE = 'stray', M, LU, LV, discard)
}

rInvWishRevCholesky_test <- function(v, Psi) {
    .Call('_stray_rInvWishRevCholesky_test', PACKAGE = 'stray', v, Psi)
}

rInvWishRevCholesky_thread_test <- function(v, Psi, discard) {
    .Call('_stray_rInvWishRevCholesky_thread_test', PACKAGE = 'stray', v, Psi, discard)
}

rInvWishRevCholesky_thread_inplace_test <- function(v, Psi, discard) {
    .Call('_stray_rInvWishRevCholesky_thread_inplace_test', PACKAGE = 'stray', v, Psi, discard)
}

rMatUnitNormal_test1 <- function(n, m) {
    .Call('_stray_rMatUnitNormal_test1', PACKAGE = 'stray', n, m)
}

rMatUnitNormal_test2 <- function(n) {
    .Call('_stray_rMatUnitNormal_test2', PACKAGE = 'stray', n)
}

#' Log of Multivarate Gamma Function - Gamma_p(a)
#' @param a defined by Gamma_p(a)
#' @param p defined by Gamma_p(a)
#' @references https://en.wikipedia.org/wiki/Multivariate_gamma_function
lmvgamma <- function(a, p) {
    .Call('_stray_lmvgamma', PACKAGE = 'stray', a, p)
}

#' Derivative of Log of Multivariate Gamma Function - Gamma_p(a)
#' @param a defined by Gamma_p(a)
#' @param p defined by Gamma_p(a)
#' @references https://en.wikipedia.org/wiki/Multivariate_gamma_function
lmvgamma_deriv <- function(a, p) {
    .Call('_stray_lmvgamma_deriv', PACKAGE = 'stray', a, p)
}

eigen_lap_test <- function(n_samples, m, S, eigvalthresh) {
    .Call('_stray_eigen_lap_test', PACKAGE = 'stray', n_samples, m, S, eigvalthresh)
}

cholesky_lap_test <- function(n_samples, m, S, eigvalthresh) {
    .Call('_stray_cholesky_lap_test', PACKAGE = 'stray', n_samples, m, S, eigvalthresh)
}

LaplaceApproximation_test <- function(n_samples, m, S, decomp_method, eigvalthresh) {
    .Call('_stray_LaplaceApproximation_test', PACKAGE = 'stray', n_samples, m, S, decomp_method, eigvalthresh)
}

alrInv_default_test <- function(eta) {
    .Call('_stray_alrInv_default_test', PACKAGE = 'stray', eta)
}

alr_default_test <- function(pi) {
    .Call('_stray_alr_default_test', PACKAGE = 'stray', pi)
}

rDirichlet_test <- function(n_samples, alpha) {
    .Call('_stray_rDirichlet_test', PACKAGE = 'stray', n_samples, alpha)
}

MultDirichletBoot_test <- function(n_samples, eta, Y, pseudocount) {
    .Call('_stray_MultDirichletBoot_test', PACKAGE = 'stray', n_samples, eta, Y, pseudocount)
}

fillUnitNormal_test <- function(Z) {
    invisible(.Call('_stray_fillUnitNormal_test', PACKAGE = 'stray', Z))
}