R/RcppExports.R

In GPBayes: Tools for Gaussian Process Modeling in Uncertainty Quantification

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

#' @title Confluent hypergeometric function of the second kind
#' @description This function calls the GSL scientific library to evaluate 
#' the confluent hypergeometric function of the second kind; see Abramowitz and Stegun 1972, p.505.
#' 
#' @param a a real value
#' @param b a real value
#' @param x a real value
#'
#'
#' @author Pulong Ma \email{mpulong@@gmail.com}
#' @seealso \code{\link{CH}}
#' @return a numerical value
HypergU <- function(a, b, x) {
    .Call(`_GPBayes_HypergU`, a, b, x)
}

#' @title Modified Bessel function of the second kind
#' @description This function calls the GSL scientific library to evaluate 
#' the modified Bessel function of the second kind.
#' 
#' @param nu a real positive value
#' @param z a real positive value
#'
#' @author Pulong Ma \email{mpulong@@gmail.com}
#' @seealso \code{\link{matern}}
#' @return a numerical value 
BesselK <- function(nu, z) {
    .Call(`_GPBayes_BesselK`, nu, z)
}

#' @title The Confluent Hypergeometric correlation function proposed by Ma and Bhadra (2023)
#'
#' @description This function computes the Confluent Hypergeometric correlation function given
#' a distance matrix. The Confluent Hypergeometric correlation function is given by 
#' \deqn{C(h) = \frac{\Gamma(\nu+\alpha)}{\Gamma(\nu)} 
#' \mathcal{U}\left(\alpha, 1-\nu,  \biggr(\frac{h}{\beta}\biggr)^2 \right),}
#' where \eqn{\alpha} is the tail decay parameter. \eqn{\beta} is the range parameter.
#' \eqn{\nu} is the smoothness parameter. \eqn{\mathcal{U}(\cdot)} is the confluent hypergeometric
#' function of the second kind. Note that this parameterization of the CH covariance
#' is different from the one in Ma and Bhadra (2023). For details about this covariance, 
#' see Ma and Bhadra (2023; \doi{10.1080/01621459.2022.2027775}).  
#' 
#' @param d a matrix of distances
#' @param range a numerical value containing the range parameter 
#' @param tail a numerical value containing the tail decay parameter
#' @param nu a numerical value containing the smoothness parameter
#'
#' @author Pulong Ma \email{mpulong@@gmail.com}
#' @return a numerical matrix  
#' @seealso \link{GPBayes-package}, \code{\link{GaSP}}, \linkS4class{gp}, \code{\link{matern}}, \code{\link{kernel}}, \code{\link{ikernel}}
CH <- function(d, range, tail, nu) {
    .Call(`_GPBayes_CH`, d, range, tail, nu)
}

CH_deriv_range <- function(d, range, tail, nu) {
    .Call(`_GPBayes_CH_deriv_range`, d, range, tail, nu)
}

CH_deriv_tail <- function(d, range, tail, nu) {
    .Call(`_GPBayes_CH_deriv_tail`, d, range, tail, nu)
}

CH_deriv_nu <- function(d, range, tail, nu) {
    .Call(`_GPBayes_CH_deriv_nu`, d, range, tail, nu)
}

deriv_ARD_CH <- function(d, range, tail, nu) {
    .Call(`_GPBayes_deriv_ARD_CH`, d, range, tail, nu)
}

#' @title The Matérn correlation function proposed by Matérn (1960)
#' 
#' @description This function computes the Matérn correlation function given
#' a distance matrix. The Matérn correlation function is given by
#' \deqn{C(h)=\frac{2^{1-\nu}}{\Gamma(\nu)} \left(\frac{h}{\phi} \right)^{\nu} 
#' \mathcal{K}_{\nu}\left( \frac{h}{\phi} \right),}
#' where \eqn{\phi} is the range parameter. \eqn{\nu} is the smoothness parameter. 
#' \eqn{\mathcal{K}_{\nu}(\cdot)} is the modified Bessel function of the second kind of order \eqn{\nu}.
#' The form of covariance includes the following special cases by specifying \eqn{\nu} to be 0.5, 1.5, 2.5.
#' \itemize{
#' \item{\eqn{\nu=0.5} corresponds to the  exponential correlation function (\strong{exp}) 
#' of the form
#' \deqn{C(h) = \exp\left\{ - \frac{h}{\phi} \right\} }
#'}
#' \item{\eqn{\nu=1.5} corresponds to the Matérn correlation function with smoothness parameter 1.5 (\strong{matern_3_2}) 
#' of the form
#' \deqn{C(h) = \left( 1 + \frac{h}{\phi} \right) \exp\left\{ - \frac{h}{\phi} \right\} }
#'}
#' \item{\eqn{\nu=2.5} corresponds to the Matérn correlation function with smoothness parameter 2.5 (\strong{matern_5_2}) 
#' of the form
#' \deqn{C(h) = \left\{ 1 + \frac{h}{\phi}  + \frac{1}{3}\left(\frac{h}{\phi}\right)^2 \right\} \exp\left\{ - \frac{h}{\phi} \right\} }
#'}
#'}
#' @param d a matrix of distances
#' @param range a numerical value containing the range parameter 
#' @param nu a numerical value containing the smoothness parameter
#'
#' @author Pulong Ma \email{mpulong@@gmail.com}
#' @return a numerical matrix 
#' @seealso \link{GPBayes-package}, \code{\link{GaSP}}, \linkS4class{gp}, \code{\link{CH}}, \code{\link{kernel}}, \code{\link{ikernel}}
#'
matern <- function(d, range, nu) {
    .Call(`_GPBayes_matern`, d, range, nu)
}

matern_deriv_range <- function(d, range, nu) {
    .Call(`_GPBayes_matern_deriv_range`, d, range, nu)
}

deriv_ARD_matern <- function(d, range, nu) {
    .Call(`_GPBayes_deriv_ARD_matern`, d, range, nu)
}

#' @title The powered-exponential correlation function
#'
#' @description This function computes the powered-exponential correlation function given
#' a distance matrix. The powered-exponential correlation function is given by
#'                \deqn{C(h)=\exp\left\{-\left(\frac{h}{\phi}\right)^{\nu}\right\},}
#' where \eqn{\phi} is the range parameter. \eqn{\nu} is the smoothness parameter.
#' The case \eqn{\nu=2} corresponds to the well-known Gaussian correlation.
#' @param d a matrix of distances
#' @param range a numerical value containing the range parameter 
#' @param nu a numerical value containing the smoothness parameter
#' 
#' @author Pulong Ma \email{mpulong@@gmail.com}
#' @return a numerical matrix 
#' @seealso \code{\link{kernel}}
#'
powexp <- function(d, range, nu) {
    .Call(`_GPBayes_powexp`, d, range, nu)
}

deriv_ARD_powexp <- function(d, range, nu) {
    .Call(`_GPBayes_deriv_ARD_powexp`, d, range, nu)
}

#' @title The generalized Cauchy correlation function
#'
#' @description This function computes the generalized Cauchy correlation function given
#' a distance matrix. The generalized Cauchy covariance is given by
#' \deqn{C(h) = \left\{ 1 + \left( \frac{h}{\phi} \right)^{\nu}  
#'             \right\}^{-\alpha/\nu},}
#' where \eqn{\phi} is the range parameter. \eqn{\alpha} is the tail decay parameter.
#' \eqn{\nu} is the smoothness parameter. 
#' The case where \eqn{\nu=2} corresponds to the Cauchy covariance model, which is infinitely differentiable.
#' 
#' @param d a matrix of distances
#' @param range a numerical value containing the range parameter 
#' @param tail a numerical value containing the tail decay parameter
#' @param nu a numerical value containing the smoothness parameter
#'
#' @author Pulong Ma \email{mpulong@@gmail.com}
#' @return a numerical matrix 
#' @seealso \code{\link{kernel}}
#'
cauchy <- function(d, range, tail, nu) {
    .Call(`_GPBayes_cauchy`, d, range, tail, nu)
}

#' @title A wraper to build different kinds of correlation matrices between two sets of inputs
#' @name ikernel
#' @description This function wraps existing built-in routines to construct a covariance 
#' matrix for two input matrices based on data type, covariance type, and distance type. The constructed 
#' covariance matrix can be directly used for GaSP fitting and and prediction for spatial 
#' data, spatio-temporal data, and computer experiments. This function explicitly takes inputs as arguments. 
#' The prefix ``i" in \code{\link{ikernel}} standards for ``input".
#'
#' @param input1 a matrix of input locations 
#' @param input2 a matrix of input locations
#' 
#' @param range a vector of range parameters, which could be a scalar.
#' @param tail a vector of tail decay parameters, which could be a scalar.
#' @param nu a vector of smoothness parameters, which could be a scalar.
#' @param covmodel a list of two strings: \strong{family}, \strong{form}, where \strong{family} indicates the family of covariance functions 
#' including the Confluent Hypergeometric class, the Matérn class, the Cauchy class, the powered-exponential class. \strong{form} indicates the 
#' specific form of covariance structures including the isotropic form, tensor form, automatic relevance determination form. 
#' \describe{
#' \item{\strong{family}}{
#' \describe{
#' \item{CH}{The Confluent Hypergeometric correlation function is given by 
#' \deqn{C(h) = \frac{\Gamma(\nu+\alpha)}{\Gamma(\nu)} 
#' \mathcal{U}\left(\alpha, 1-\nu, \left(\frac{h}{\beta}\right)^2\right),}
#' where \eqn{\alpha} is the tail decay parameter. \eqn{\beta} is the range parameter.
#' \eqn{\nu} is the smoothness parameter. \eqn{\mathcal{U}(\cdot)} is the confluent hypergeometric
#' function of the second kind. For details about this covariance, 
#' see Ma and Bhadra (2023; \doi{10.1080/01621459.2022.2027775}).  
#' }
#' \item{cauchy}{The generalized Cauchy covariance is given by
#' \deqn{C(h) = \left\{ 1 + \left( \frac{h}{\phi} \right)^{\nu}  
#'             \right\}^{-\alpha/\nu},}
#' where \eqn{\phi} is the range parameter. \eqn{\alpha} is the tail decay parameter.
#' \eqn{\nu} is the smoothness parameter with default value at 2.
#'}
#'
#' \item{matern}{The Matérn correlation function is given by
#' \deqn{C(h)=\frac{2^{1-\nu}}{\Gamma(\nu)} \left( \frac{h}{\phi} \right)^{\nu} 
#' \mathcal{K}_{\nu}\left( \frac{h}{\phi} \right),}
#' where \eqn{\phi} is the range parameter. \eqn{\nu} is the smoothness parameter. 
#' \eqn{\mathcal{K}_{\nu}(\cdot)} is the modified Bessel function of the second kind of order \eqn{\nu}.
#' }
#' \describe{
#' \item{exp}{This is the Matérn correlation with \eqn{\nu=0.5}. This covariance should be specified as \strong{matern} with smoothness parameter \eqn{\nu=0.5}.
#' }
#' \item{matern_3_2}{This is the Matérn correlation with \eqn{\nu=1.5}.
#' This covariance should be specified as \strong{matern} with smoothness parameter \eqn{\nu=1.5}.}
#' \item{matern_5_2}{This is the Matérn correlation with \eqn{\nu=2.5}. 
#' This covariance should be specified as \strong{matern} with smoothness parameter \eqn{\nu=2.5}.}
#' }
#'
#'
#' \item{powexp}{The powered-exponential correlation function is given by
#'                \deqn{C(h)=\exp\left\{-\left(\frac{h}{\phi}\right)^{\nu}\right\},}
#' where \eqn{\phi} is the range parameter. \eqn{\nu} is the smoothness parameter.
#' }
#' \item{gauss}{The Gaussian correlation function is given by 
#' \deqn{C(h)=\exp\left(-\frac{h^2}{\phi^2}\right),}
#' where \eqn{\phi} is the range parameter.
#'  }
#' }
#' }
#' 
#' \item{\strong{form}}{
#' \describe{
#'  \item{isotropic}{This indicates the isotropic form of covariance functions. That is,
#'  \deqn{C(\mathbf{h}) = C^0(\|\mathbf{h}\|; \boldsymbol \theta),} where \eqn{\| \mathbf{h}\|} denotes the 
#' Euclidean distance or the great circle distance for data on sphere. \eqn{C^0(\cdot)} denotes 
#' any isotropic covariance family specified in \strong{family}.}
#'  \item{tensor}{This indicates the tensor product of correlation functions. That is, 
#' \deqn{ C(\mathbf{h}) = \prod_{i=1}^d C^0(|h_i|; \boldsymbol \theta_i),}
#' where \eqn{d} is the dimension of input space. \eqn{h_i} is the distance along the \eqn{i}th input dimension. This type of covariance structure has been often used in Gaussian process emulation for computer experiments.
#'}
#'  \item{ARD}{This indicates the automatic relevance determination form. That is, 
#' \deqn{C(\mathbf{h}) = C^0\left(\sqrt{\sum_{i=1}^d\frac{h_i^2}{\phi^2_i}}; \boldsymbol \theta \right),}
#' where \eqn{\phi_i} denotes the range parameter along the \eqn{i}th input dimension.}
#'  }
#' }
#'
#'}
#' 
#' @param dtype a string indicating distance type: \strong{Euclidean}, \strong{GCD}, where the latter indicates great circle distance.
NULL

#' @title A wraper to compute the natural logarithm of the integrated likelihood function
#' 
#' @description This function wraps existing built-in routines to construct  
#' the natural logarithm of the integrated likelihood function. The constructed 
#' loglikelihood can be directly used for numerical optimization 
#'
#' @param par a numerical vector, with which numerical optimization routine such as \code{\link[stats]{optim}} can be
#' carried out directly. When the confluent Hypergeometric class is used, it is used to hold values 
#' for \strong{range}, \strong{tail}, \strong{nugget}, and \strong{nu} if the smoothness parameter is estimated.
#' When the Matérn class or powered-exponential class is used, it is used to hold values
#' for \strong{range}, \strong{nugget}, and \strong{nu} if the smoothness parameter is estimated.
#' The order of the parameter values in \code{par} cannot be changed. For tensor or ARD form correlation
#' functions, \strong{range} and \strong{tail} becomes a vector. 
#' @param output a matrix of outputs 
#' @param H a matrix of regressors in the mean function of a GaSP model.
#' @param d an R object holding the distances. It should be a distance matrix for constructing isotropic correlation matrix,
#' or a list of distance matrices along each input dimension for constructing tensor or ARD types of correlation matrix.
#' @param covmodel a list of two strings: \strong{family}, \strong{form}, where \strong{family} indicates the family of covariance functions 
#' including the Confluent Hypergeometric class, the Matérn class, the Cauchy class, the powered-exponential class. \strong{form} indicates the 
#' specific form of covariance structures including the isotropic form, tensor form, automatic relevance determination form. 
#' \describe{
#' \item{\strong{family}}{
#' \describe{
#' \item{CH}{The Confluent Hypergeometric correlation function is given by 
#' \deqn{C(h) = \frac{\Gamma(\nu+\alpha)}{\Gamma(\nu)} 
#' \mathcal{U}\left(\alpha, 1-\nu, \left(\frac{h}{\beta}\right)^2\right),}
#' where \eqn{\alpha} is the tail decay parameter. \eqn{\beta} is the range parameter.
#' \eqn{\nu} is the smoothness parameter. \eqn{\mathcal{U}(\cdot)} is the confluent hypergeometric
#' function of the second kind. For details about this covariance, 
#' see Ma and Bhadra (2023; \doi{10.1080/01621459.2022.2027775}).  
#' }
#' \item{cauchy}{The generalized Cauchy covariance is given by
#' \deqn{C(h) = \left\{ 1 + \left( \frac{h}{\phi} \right)^{\nu}  
#'             \right\}^{-\alpha/\nu},}
#' where \eqn{\phi} is the range parameter. \eqn{\alpha} is the tail decay parameter.
#' \eqn{\nu} is the smoothness parameter with default value at 2.
#'}
#'
#' \item{matern}{The Matérn correlation function is given by
#' \deqn{C(h)=\frac{2^{1-\nu}}{\Gamma(\nu)} \left( \frac{h}{\phi} \right)^{\nu} 
#' \mathcal{K}_{\nu}\left( \frac{h}{\phi} \right),}
#' where \eqn{\phi} is the range parameter. \eqn{\nu} is the smoothness parameter. 
#' \eqn{\mathcal{K}_{\nu}(\cdot)} is the modified Bessel function of the second kind of order \eqn{\nu}.
#' }
#' \describe{
#' \item{exp}{This is the Matérn correlation with \eqn{\nu=0.5}. This covariance should be specified as \strong{matern} with smoothness parameter \eqn{\nu=0.5}.
#' }
#' \item{matern_3_2}{This is the Matérn correlation with \eqn{\nu=1.5}.
#' This covariance should be specified as \strong{matern} with smoothness parameter \eqn{\nu=1.5}.}
#' \item{matern_5_2}{This is the Matérn correlation with \eqn{\nu=2.5}. 
#' This covariance should be specified as \strong{matern} with smoothness parameter \eqn{\nu=2.5}.}
#' }
#'
#'
#' \item{powexp}{The powered-exponential correlation function is given by
#'                \deqn{C(h)=\exp\left\{-\left(\frac{h}{\phi}\right)^{\nu}\right\},}
#' where \eqn{\phi} is the range parameter. \eqn{\nu} is the smoothness parameter.
#' }
#' \item{gauss}{The Gaussian correlation function is given by 
#' \deqn{C(h)=\exp\left(-\frac{h^2}{\phi^2}\right),}
#' where \eqn{\phi} is the range parameter.
#'  }
#' }
#' }
#' 
#' \item{\strong{form}}{
#' \describe{
#'  \item{isotropic}{This indicates the isotropic form of covariance functions. That is,
#'  \deqn{C(\mathbf{h}) = C^0(\|\mathbf{h}\|; \boldsymbol \theta),} where \eqn{\| \mathbf{h}\|} denotes the 
#' Euclidean distance or the great circle distance for data on sphere. \eqn{C^0(\cdot)} denotes 
#' any isotropic covariance family specified in \strong{family}.}
#'  \item{tensor}{This indicates the tensor product of correlation functions. That is, 
#' \deqn{ C(\mathbf{h}) = \prod_{i=1}^d C^0(|h_i|; \boldsymbol \theta_i),}
#' where \eqn{d} is the dimension of input space. \eqn{h_i} is the distance along the \eqn{i}th input dimension. This type of covariance structure has been often used in Gaussian process emulation for computer experiments.
#'}
#'  \item{ARD}{This indicates the automatic relevance determination form. That is, 
#' \deqn{C(\mathbf{h}) = C^0\left(\sqrt{\sum_{i=1}^d\frac{h_i^2}{\phi^2_i}}; \boldsymbol \theta \right),}
#' where \eqn{\phi_i} denotes the range parameter along the \eqn{i}th input dimension.}
#'  }
#' }
#'
#'}
#' 
#' @param smooth The smoothness parameter \eqn{\nu} in a correlation function.
#' @param smoothness_est a logical value indicating whether the smoothness parameter is estimated. 
#' @return The natural logarithm of marginal or integrated likelihood
#' @author Pulong Ma \email{mpulong@@gmail.com}
#' @seealso \code{\link{CH}}, \code{\link{matern}}, \code{\link{gp.optim}}, \link{GPBayes-package}, \code{\link{GaSP}}
loglik <- function(par, output, H, d, covmodel, smooth, smoothness_est) {
    .Call(`_GPBayes_loglik`, par, output, H, d, covmodel, smooth, smoothness_est)
}

iso_loglik <- function(par, output, H, d, covmodel, smooth, smoothness_est) {
    .Call(`_GPBayes_iso_loglik`, par, output, H, d, covmodel, smooth, smoothness_est)
}

tensor_loglik <- function(par, output, H, d, covmodel, smooth, smoothness_est) {
    .Call(`_GPBayes_tensor_loglik`, par, output, H, d, covmodel, smooth, smoothness_est)
}

ARD_loglik <- function(par, output, H, d, covmodel, smooth, smoothness_est) {
    .Call(`_GPBayes_ARD_loglik`, par, output, H, d, covmodel, smooth, smoothness_est)
}

gradient_loglik <- function(par, output, H, d, covmodel, smooth, smoothness_est) {
    .Call(`_GPBayes_gradient_loglik`, par, output, H, d, covmodel, smooth, smoothness_est)
}

loglik_xi <- function(par, output, H, d, covmodel, smooth, smoothness_est) {
    .Call(`_GPBayes_loglik_xi`, par, output, H, d, covmodel, smooth, smoothness_est)
}

gradient_loglik_xi <- function(par, output, H, d, covmodel, smooth, smoothness_est) {
    .Call(`_GPBayes_gradient_loglik_xi`, par, output, H, d, covmodel, smooth, smoothness_est)
}

#' @title Compute distances for two sets of inputs 
#' 
#' @description This function computes distances for two sets of inputs and returns
#' a \code{R} object. 
#'
#' @param input1 a matrix of inputs
#' @param input2 a matrix of inputs
#' @param type a string indicating the form of distances with three froms supported currently: \strong{isotropic}, \strong{tensor}, \strong{ARD}.
#' @param dtype a string indicating distance type: \strong{Euclidean}, \strong{GCD}, where the latter indicates great circle distance.
#' @return a R object holding distances for two sets of inputs. If \strong{type} is \strong{isotropic}, a matrix of distances
#' is returned; if \strong{type} is \strong{tensor} or \strong{ARD}, a list of distance matrices 
#' along each input dimension is returned.
#' @author Pulong Ma \email{mpulong@@gmail.com}
#' @return a numeric vector or matrix of distances 
#' @examples
#' input = seq(0,1,length=20)
#' d = distance(input, input, type="isotropic", dtype="Euclidean")
#'
distance <- function(input1, input2, type = "isotropic", dtype = "Euclidean") {
    .Call(`_GPBayes_distance`, input1, input2, type, dtype)
}

#' @return a correlation matrix
#' @author Pulong Ma \email{mpulong@@gmail.com}
#' @examples
#' input = seq(0,1,length=10)
#' 
#' cormat = ikernel(input,input,range=0.5,tail=0.2,nu=2.5,
#'          covmodel=list(family="CH",form="isotropic"))
#'
#' @seealso \code{\link{CH}}, \code{\link{matern}}, \code{\link{kernel}}, \link{GPBayes-package}, \code{\link{GaSP}}
ikernel <- function(input1, input2, range, tail, nu, covmodel, dtype = "Euclidean") {
    .Call(`_GPBayes_ikernel`, input1, input2, range, tail, nu, covmodel, dtype)
}

#' @title A wraper to build different kinds of correlation matrices with distance as arguments 
#' 
#' @description This function wraps existing built-in routines to construct a covariance 
#' matrix based on data type, covariance type, and distance type with distances as inputs. The constructed 
#' covariance matrix can be directly used for GaSP fitting and and prediction for spatial 
#' data, spatio-temporal data, and computer experiments. 
#'
#' @param d a matrix or a list of distances
#' 
#' @param range a vector of range parameters, which could be a scalar. 
#' @param tail a vector of tail decay parameters, which could be a scalar.
#' @param nu a vector of smoothness parameters, which could be a scalar.
#' @param covmodel a list of two strings: \strong{family}, \strong{form}, where \strong{family} indicates the family of covariance functions 
#' including the Confluent Hypergeometric class, the Matérn class, the Cauchy class, the powered-exponential class. \strong{form} indicates the 
#' specific form of covariance structures including the isotropic form, tensor form, automatic relevance determination form. 
#' \describe{
#' \item{\strong{family}}{
#' \describe{
#' \item{CH}{The Confluent Hypergeometric correlation function is given by 
#' \deqn{C(h) = \frac{\Gamma(\nu+\alpha)}{\Gamma(\nu)} 
#' \mathcal{U}\left(\alpha, 1-\nu, \left(\frac{h}{\beta}\right)^2\right),}
#' where \eqn{\alpha} is the tail decay parameter. \eqn{\beta} is the range parameter.
#' \eqn{\nu} is the smoothness parameter. \eqn{\mathcal{U}(\cdot)} is the confluent hypergeometric
#' function of the second kind. For details about this covariance, 
#' see Ma and Bhadra (2023; \doi{10.1080/01621459.2022.2027775}).  
#' }
#' \item{cauchy}{The generalized Cauchy covariance is given by
#' \deqn{C(h) = \left\{ 1 + \left( \frac{h}{\phi} \right)^{\nu}  
#'             \right\}^{-\alpha/\nu},}
#' where \eqn{\phi} is the range parameter. \eqn{\alpha} is the tail decay parameter.
#' \eqn{\nu} is the smoothness parameter with default value at 2.
#'}
#'
#' \item{matern}{The Matérn correlation function is given by
#' \deqn{C(h)=\frac{2^{1-\nu}}{\Gamma(\nu)} \left( \frac{h}{\phi} \right)^{\nu} 
#' \mathcal{K}_{\nu}\left( \frac{h}{\phi} \right),}
#' where \eqn{\phi} is the range parameter. \eqn{\nu} is the smoothness parameter. 
#' \eqn{\mathcal{K}_{\nu}(\cdot)} is the modified Bessel function of the second kind of order \eqn{\nu}.
#' }
#' \describe{
#' \item{exp}{This is the Matérn correlation with \eqn{\nu=0.5}. This covariance should be specified as \strong{matern} with smoothness parameter \eqn{\nu=0.5}.
#' }
#' \item{matern_3_2}{This is the Matérn correlation with \eqn{\nu=1.5}.
#' This covariance should be specified as \strong{matern} with smoothness parameter \eqn{\nu=1.5}.}
#' \item{matern_5_2}{This is the Matérn correlation with \eqn{\nu=2.5}. 
#' This covariance should be specified as \strong{matern} with smoothness parameter \eqn{\nu=2.5}.}
#' }
#'
#' \item{powexp}{The powered-exponential correlation function is given by
#'                \deqn{C(h)=\exp\left\{-\left(\frac{h}{\phi}\right)^{\nu}\right\},}
#' where \eqn{\phi} is the range parameter. \eqn{\nu} is the smoothness parameter.
#' }
#' \item{gauss}{The Gaussian correlation function is given by 
#' \deqn{C(h)=\exp\left(-\frac{h^2}{\phi^2}\right),}
#' where \eqn{\phi} is the range parameter.
#'  }
#' }
#' }
#' 
#' \item{\strong{form}}{
#' \describe{
#'  \item{isotropic}{This indicates the isotropic form of covariance functions. That is,
#'  \deqn{C(\mathbf{h}) = C^0(\|\mathbf{h}\|; \boldsymbol \theta),} where \eqn{\| \mathbf{h}\|} denotes the 
#' Euclidean distance or the great circle distance for data on sphere. \eqn{C^0(\cdot)} denotes 
#' any isotropic covariance family specified in \strong{family}.}
#'  \item{tensor}{This indicates the tensor product of correlation functions. That is, 
#' \deqn{ C(\mathbf{h}) = \prod_{i=1}^d C^0(|h_i|; \boldsymbol \theta_i),}
#' where \eqn{d} is the dimension of input space. \eqn{h_i} is the distance along the \eqn{i}th input dimension. This type of covariance structure has been often used in Gaussian process emulation for computer experiments.
#'}
#'  \item{ARD}{This indicates the automatic relevance determination form. That is, 
#' \deqn{C(\mathbf{h}) = C^0\left(\sqrt{\sum_{i=1}^d\frac{h_i^2}{\phi^2_i}}; \boldsymbol \theta \right),}
#' where \eqn{\phi_i} denotes the range parameter along the \eqn{i}th input dimension.}
#'  }
#' }
#'
#'}
#' 
#' @return a correlation matrix 
#' @author Pulong Ma \email{mpulong@@gmail.com}
#' @seealso \code{\link{CH}}, \code{\link{matern}}, \code{\link{ikernel}}, \link{GPBayes-package}, \code{\link{GaSP}}
#' @examples
#' input = seq(0,1,length=10)
#' d = distance(input,input,type="isotropic",dtype="Euclidean")
#' cormat = kernel(d,range=0.5,tail=0.2,nu=2.5,
#'          covmodel=list(family="CH",form="isotropic"))
kernel <- function(d, range, tail, nu, covmodel) {
    .Call(`_GPBayes_kernel`, d, range, tail, nu, covmodel)
}

#' @title A wraper to construct the derivative of correlation matrix with respect to correlation parameters 
#' 
#' @description This function wraps existing built-in routines to construct the 
#' derivative of correlation matrix with respect to correlation parameters. 
#' @param d a matrix or a list of distances returned from \code{\link{distance}}.
#' 
#' @param range a vector of range parameters 
#' @param tail a vector of tail decay parameters
#' @param nu a vector of smoothness parameters
#' @param covmodel a list of two strings: \strong{family}, \strong{form}, where \strong{family} indicates the family of covariance functions 
#' including the Confluent Hypergeometric class, the Matérn class, the Cauchy class, the powered-exponential class. \strong{form} indicates the 
#' specific form of covariance structures including the isotropic form, tensor form, automatic relevance determination form. 
#' \describe{
#' \item{\strong{family}}{
#' \describe{
#' \item{CH}{The Confluent Hypergeometric correlation function is given by 
#' \deqn{C(h) = \frac{\Gamma(\nu+\alpha)}{\Gamma(\nu)} 
#' \mathcal{U}\left(\alpha, 1-\nu, \left(\frac{h}{\beta}\right)^2\right),}
#' where \eqn{\alpha} is the tail decay parameter. \eqn{\beta} is the range parameter.
#' \eqn{\nu} is the smoothness parameter. \eqn{\mathcal{U}(\cdot)} is the confluent hypergeometric
#' function of the second kind. For details about this covariance, 
#' see Ma and Bhadra (2023; \doi{10.1080/01621459.2022.2027775}).  
#' }
#' \item{cauchy}{The generalized Cauchy covariance is given by
#' \deqn{C(h) = \left\{ 1 + \left( \frac{h}{\phi} \right)^{\nu}  
#'             \right\}^{-\alpha/\nu},}
#' where \eqn{\phi} is the range parameter. \eqn{\alpha} is the tail decay parameter.
#' \eqn{\nu} is the smoothness parameter with default value at 2.
#'}
#'
#' \item{matern}{The Matérn correlation function is given by
#' \deqn{C(h)=\frac{2^{1-\nu}}{\Gamma(\nu)} \left( \frac{h}{\phi} \right)^{\nu} 
#' \mathcal{K}_{\nu}\left( \frac{h}{\phi} \right),}
#' where \eqn{\phi} is the range parameter. \eqn{\nu} is the smoothness parameter. 
#' \eqn{\mathcal{K}_{\nu}(\cdot)} is the modified Bessel function of the second kind of order \eqn{\nu}.
#' }
#' \describe{
#' \item{exp}{This is the Matérn correlation with \eqn{\nu=0.5}. This covariance should be specified as \strong{matern} with smoothness parameter \eqn{\nu=0.5}.
#' }
#' \item{matern_3_2}{This is the Matérn correlation with \eqn{\nu=1.5}.
#' This covariance should be specified as \strong{matern} with smoothness parameter \eqn{\nu=1.5}.}
#' \item{matern_5_2}{This is the Matérn correlation with \eqn{\nu=2.5}. 
#' This covariance should be specified as \strong{matern} with smoothness parameter \eqn{\nu=2.5}.}
#' }
#'
#'
#' \item{powexp}{The powered-exponential correlation function is given by
#'                \deqn{C(h)=\exp\left\{-\left(\frac{h}{\phi}\right)^{\nu}\right\},}
#' where \eqn{\phi} is the range parameter. \eqn{\nu} is the smoothness parameter.
#' }
#' \item{gauss}{The Gaussian correlation function is given by 
#' \deqn{C(h)=\exp\left(-\frac{h^2}{\phi^2}\right),}
#' where \eqn{\phi} is the range parameter.
#'  }
#' }
#' }
#' 
#' \item{\strong{form}}{
#' \describe{
#'  \item{isotropic}{This indicates the isotropic form of covariance functions. That is,
#'  \deqn{C(\mathbf{h}) = C^0(\|\mathbf{h}\|; \boldsymbol \theta),} where \eqn{\| \mathbf{h}\|} denotes the 
#' Euclidean distance or the great circle distance for data on sphere. \eqn{C^0(\cdot)} denotes 
#' any isotropic covariance family specified in \strong{family}.}
#'  \item{tensor}{This indicates the tensor product of correlation functions. That is, 
#' \deqn{ C(\mathbf{h}) = \prod_{i=1}^d C^0(|h_i|; \boldsymbol \theta_i),}
#' where \eqn{d} is the dimension of input space. \eqn{h_i} is the distance along the \eqn{i}th input dimension. This type of covariance structure has been often used in Gaussian process emulation for computer experiments.
#'}
#'  \item{ARD}{This indicates the automatic relevance determination form. That is, 
#' \deqn{C(\mathbf{h}) = C^0\left(\sqrt{\sum_{i=1}^d\frac{h_i^2}{\phi^2_i}}; \boldsymbol \theta \right),}
#' where \eqn{\phi_i} denotes the range parameter along the \eqn{i}th input dimension.}
#'  }
#' }
#'
#'}
#' 
#' @return a list of matrices 
#' @author Pulong Ma \email{mpulong@@gmail.com}
#' @seealso \code{\link{CH}}, \code{\link{matern}}, \code{\link{kernel}}, \link{GPBayes-package}, \code{\link{GaSP}}
#' @examples
#' input = seq(0,1,length=10)
#' d = distance(input,input,type="isotropic",dtype="Euclidean")
#' dR = deriv_kernel(d,range=0.5,tail=0.2,nu=2.5,
#'          covmodel=list(family="CH",form="isotropic"))
deriv_kernel <- function(d, range, tail, nu, covmodel) {
    .Call(`_GPBayes_deriv_kernel`, d, range, tail, nu, covmodel)
}

reference_prior <- function(H, d, par, covmodel, smoothness_est) {
    .Call(`_GPBayes_reference_prior`, H, d, par, covmodel, smoothness_est)
}

FisherInfo <- function(input, sig2, range, tail, nu, nugget, covmodel, dtype) {
    .Call(`_GPBayes_FisherInfo`, input, sig2, range, tail, nu, nugget, covmodel, dtype)
}

FisherIR_intlik <- function(H, input, range, tail, nu, nugget, covmodel, dtype) {
    .Call(`_GPBayes_FisherIR_intlik`, H, input, range, tail, nu, nugget, covmodel, dtype)
}

GPsim <- function(input, H, par, covmodel, nsample = 1L, dtype = "Euclidean") {
    .Call(`_GPBayes_GPsim`, input, H, par, covmodel, nsample, dtype)
}

GPpredict <- function(output, H, input, input_new, Hnew, par, covmodel, dtype) {
    .Call(`_GPBayes_GPpredict`, output, H, input, input_new, Hnew, par, covmodel, dtype)
}

post_predictive_sampling <- function(output, H, input, input_new, Hnew, MCMCsample, smooth, smoothness_est, covmodel, dtype) {
    .Call(`_GPBayes_post_predictive_sampling`, output, H, input, input_new, Hnew, MCMCsample, smooth, smoothness_est, covmodel, dtype)
}

condsim <- function(output, H, input, input_new, Hnew, par, covmodel, dtype, nsample) {
    .Call(`_GPBayes_condsim`, output, H, input, input_new, Hnew, par, covmodel, dtype, nsample)
}

MCMCOBayes <- function(output, H, input, par, covmodel, smoothness_est, proposal, nsample, dtype, verbose = TRUE) {
    .Call(`_GPBayes_MCMCOBayes`, output, H, input, par, covmodel, smoothness_est, proposal, nsample, dtype, verbose)
}

MCMCOBayes_pred <- function(output, H, input, input_new, Hnew, par, covmodel, smoothness_est, proposal, nsample, dtype, verbose) {
    .Call(`_GPBayes_MCMCOBayes_pred`, output, H, input, input_new, Hnew, par, covmodel, smoothness_est, proposal, nsample, dtype, verbose)
}

MCMCOBayesRef <- function(output, H, input, par, covmodel, smoothness_est, proposal, nsample, dtype, verbose) {
    .Call(`_GPBayes_MCMCOBayesRef`, output, H, input, par, covmodel, smoothness_est, proposal, nsample, dtype, verbose)
}

MCMCOBayesRef_pred <- function(output, H, input, input_new, Hnew, par, covmodel, smoothness_est, proposal, nsample, dtype, verbose) {
    .Call(`_GPBayes_MCMCOBayesRef_pred`, output, H, input, input_new, Hnew, par, covmodel, smoothness_est, proposal, nsample, dtype, verbose)
}

MCMCSBayes <- function(output, H, input, par, covmodel, smoothness_est, prior, proposal, nsample, dtype, verbose) {
    .Call(`_GPBayes_MCMCSBayes`, output, H, input, par, covmodel, smoothness_est, prior, proposal, nsample, dtype, verbose)
}

MCMCSBayes_pred <- function(output, H, input, input_new, Hnew, par, covmodel, smoothness_est, prior, proposal, nsample, dtype, verbose) {
    .Call(`_GPBayes_MCMCSBayes_pred`, output, H, input, input_new, Hnew, par, covmodel, smoothness_est, prior, proposal, nsample, dtype, verbose)
}

SPLoglik <- function(range, tail, nu, nugget, y, H, d, covmodel) {
    .Call(`_GPBayes_SPLoglik`, range, tail, nu, nugget, y, H, d, covmodel)
}

MCMCtest <- function(output, H, input, par, smoothness_est, proposal, nsample, dtype, verbose) {
    .Call(`_GPBayes_MCMCtest`, output, H, input, par, smoothness_est, proposal, nsample, dtype, verbose)
}

MCMCOBayes_Ref <- function(output, H, input, par, covmodel, smoothness_est, proposal, nsample, dtype, verbose) {
    .Call(`_GPBayes_MCMCOBayes_Ref`, output, H, input, par, covmodel, smoothness_est, proposal, nsample, dtype, verbose)
}

model_evaluation <- function(output, H, input, covmodel, smoothness, sample, output_new, input_new, Hnew, dtype, pointwise = TRUE, joint = TRUE) {
    .Call(`_GPBayes_model_evaluation`, output, H, input, covmodel, smoothness, sample, output_new, input_new, Hnew, dtype, pointwise, joint)
}