R/RcppExports.R

In GpGp: Fast Gaussian Process Computation Using Vecchia's Approximation

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

#' compute gradient of spherical harmonics functions
#' 
#' @param xyz xyz coordinates of locations on sphere
#' @param Lmax largest degree of spherical harmonics. 
#' Current only Lmax=2 supported
sph_grad_xyz <- function(xyz, Lmax) {
    .Call('_GpGp_sph_grad_xyz', PACKAGE = 'GpGp', xyz, Lmax)
}

#' Isotropic exponential covariance function
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs is a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range, nugget)
#' = \eqn{(\sigma^2,\alpha,\tau^2)}, and the covariance function is parameterized
#' as
#' \deqn{ M(x,y) = \sigma^2 exp( - || x - y ||/ \alpha )}
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
exponential_isotropic <- function(covparms, locs) {
    .Call('_GpGp_exponential_isotropic', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn exponential_isotropic Derivatives of isotropic exponential covariance
d_exponential_isotropic <- function(covparms, locs) {
    .Call('_GpGp_d_exponential_isotropic', PACKAGE = 'GpGp', covparms, locs)
}

#' Isotropic Matern covariance function
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range, smoothness, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs gives a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range, smoothness, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range, smoothness, nugget)
#' = \eqn{(\sigma^2,\alpha,\nu,\tau^2)}, and the covariance function is parameterized
#' as
#' \deqn{ M(x,y) = \sigma^2 2^{1-\nu}/\Gamma(\nu) (|| x - y ||/\alpha )^\nu K_\nu(|| x - y ||/\alpha ) }
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern_isotropic <- function(covparms, locs) {
    .Call('_GpGp_matern_isotropic', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern_isotropic Derivatives of isotropic Matern covariance
d_matern_isotropic <- function(covparms, locs) {
    .Call('_GpGp_d_matern_isotropic', PACKAGE = 'GpGp', covparms, locs)
}

#' Isotropic Matern covariance function, smoothness = 1.5
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs is a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range, nugget)
#' = \eqn{(\sigma^2,\alpha,\tau^2)}, and the covariance function is parameterized
#' as
#' \deqn{ M(x,y) = \sigma^2 (1 + || x - y || ) exp( - || x - y ||/ \alpha )}
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern15_isotropic <- function(covparms, locs) {
    .Call('_GpGp_matern15_isotropic', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn exponential_isotropic Derivatives of isotropic 
#' matern covariance with smoothness 1.5
d_matern15_isotropic <- function(covparms, locs) {
    .Call('_GpGp_d_matern15_isotropic', PACKAGE = 'GpGp', covparms, locs)
}

#' Geometrically anisotropic Matern covariance function (two dimensions)
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, L11, L21, L22, smoothness, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{2} columns.
#' Each row of locs is a point in R^2.
#' @param covparms A vector with covariance parameters
#' in the form (variance, L11, L21, L22, smoothness, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, L11, L21, L22, smoothness, nugget)
#' where L11, L21, L22, are the three non-zero entries of a lower-triangular
#' matrix L. The covariances are 
#' \deqn{ M(x,y) = \sigma^2 2^{1-\nu}/\Gamma(\nu) (|| L x - L y || )^\nu K_\nu(|| L x - L y ||) }
#' This means that L11 is interpreted as an inverse range parameter in the
#' first dimension.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern_anisotropic2D <- function(covparms, locs) {
    .Call('_GpGp_matern_anisotropic2D', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern_anisotropic2D Derivatives of anisotropic Matern covariance
d_matern_anisotropic2D <- function(covparms, locs) {
    .Call('_GpGp_d_matern_anisotropic2D', PACKAGE = 'GpGp', covparms, locs)
}

#' Geometrically anisotropic Matern covariance function (three dimensions)
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, L11, L21, L22, L31, L32, L33, smoothness, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{3} columns.
#' Each row of locs is a point in R^3.
#' @param covparms A vector with covariance parameters
#' in the form (variance, L11, L21, L22, L31, L32, L33, smoothness, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, L11, L21, L22, L31, L32, L33, smoothness, nugget)
#' where L11, L21, L22, L31, L32, L33 are the six non-zero entries of a lower-triangular
#' matrix L. The covariances are 
#' \deqn{ M(x,y) = \sigma^2 2^{1-\nu}/\Gamma(\nu) (|| L x - L y || )^\nu K_\nu(|| L x - L y ||) }
#' This means that L11 is interpreted as an inverse range parameter in the
#' first dimension.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern_anisotropic3D <- function(covparms, locs) {
    .Call('_GpGp_matern_anisotropic3D', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern_anisotropic3D Derivatives of anisotropic Matern covariance
d_matern_anisotropic3D <- function(covparms, locs) {
    .Call('_GpGp_d_matern_anisotropic3D', PACKAGE = 'GpGp', covparms, locs)
}

#' Geometrically anisotropic Matern covariance function (three dimensions, alternate parameterization)
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, B11, B12, B13, B22, B23, B33, smoothness, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{3} columns.
#' Each row of locs is a point in R^3.
#' @param covparms A vector with covariance parameters
#' in the form (variance, B11, B12, B13, B22, B23, B33, smoothness, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, B11, B12, B13, B22, B23, B33, smoothness, nugget)
#' where B11, B12, B13, B22, B23, B33, transform the three coordinates as
#' \deqn{ u_1 = B11[ x_1 + B12 x_2 + (B13 + B12 B23) x_3] }
#' \deqn{ u_2 = B22[ x_2 + B23 x_3] }
#' \deqn{ u_3 = B33[ x_3 ] }
#' NOTE: the u_1 transformation is different from previous versions of this function.
#' NOTE: now (B13,B23) can be interpreted as a drift vector in space over time.
#' Assuming x is transformed to u and y transformed to v, the covariances are 
#' \deqn{ M(x,y) = \sigma^2 2^{1-\nu}/\Gamma(\nu) (|| u - v || )^\nu K_\nu(|| u - v ||) }
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern_anisotropic3D_alt <- function(covparms, locs) {
    .Call('_GpGp_matern_anisotropic3D_alt', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern_anisotropic3D Derivatives of anisotropic Matern covariance
d_matern_anisotropic3D_alt <- function(covparms, locs) {
    .Call('_GpGp_d_matern_anisotropic3D_alt', PACKAGE = 'GpGp', covparms, locs)
}

#' Geometrically anisotropic exponential covariance function (two dimensions)
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, L11, L21, L22, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{2} columns.
#' Each row of locs is a point in R^2.
#' @param covparms A vector with covariance parameters
#' in the form (variance, L11, L21, L22, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, L11, L21, L22, nugget)
#' where L11, L21, L22, are the three non-zero entries of a lower-triangular
#' matrix L. The covariances are 
#' \deqn{ M(x,y) = \sigma^2 exp(-|| L x - L y || ) }
#' This means that L11 is interpreted as an inverse range parameter in the
#' first dimension.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
exponential_anisotropic2D <- function(covparms, locs) {
    .Call('_GpGp_exponential_anisotropic2D', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn exponential_anisotropic2D Derivatives of anisotropic exponential covariance
d_exponential_anisotropic2D <- function(covparms, locs) {
    .Call('_GpGp_d_exponential_anisotropic2D', PACKAGE = 'GpGp', covparms, locs)
}

#' Geometrically anisotropic exponential covariance function (three dimensions)
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, L11, L21, L22, L31, L32, L33, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{3} columns.
#' Each row of locs is a point in R^3.
#' @param covparms A vector with covariance parameters
#' in the form (variance, L11, L21, L22, L31, L32, L33, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, L11, L21, L22, L31, L32, L33, nugget)
#' where L11, L21, L22, L31, L32, L33 are the six non-zero entries of a lower-triangular
#' matrix L. The covariances are 
#' \deqn{ M(x,y) = \sigma^2 exp(-|| L x - L y || ) }
#' This means that L11 is interpreted as an inverse range parameter in the
#' first dimension.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
exponential_anisotropic3D <- function(covparms, locs) {
    .Call('_GpGp_exponential_anisotropic3D', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn exponential_anisotropic3D Derivatives of anisotropic exponential covariance
d_exponential_anisotropic3D <- function(covparms, locs) {
    .Call('_GpGp_d_exponential_anisotropic3D', PACKAGE = 'GpGp', covparms, locs)
}

#' Geometrically anisotropic exponential covariance function (three dimensions, alternate parameterization)
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, B11, B12, B13, B22, B23, B33, smoothness, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{3} columns.
#' Each row of locs is a point in R^3.
#' @param covparms A vector with covariance parameters
#' in the form (variance, B11, B12, B13, B22, B23, B33, smoothness, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, B11, B12, B13, B22, B23, B33, smoothness, nugget)
#' where B11, B12, B13, B22, B23, B33, transform the three coordinates as
#' \deqn{ u_1 = B11[ x_1 + B12 x_2 + (B13 + B12 B23) x_3] }
#' \deqn{ u_2 = B22[ x_2 + B23 x_3] }
#' \deqn{ u_3 = B33[ x_3 ] }
#' (B13,B23) can be interpreted as a drift vector in space over time
#' if first two dimensions are space and third is time.
#' Assuming x is transformed to u and y transformed to v, the covariances are 
#' \deqn{ M(x,y) = \sigma^2 exp( - || u - v || )
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
exponential_anisotropic3D_alt <- function(covparms, locs) {
    .Call('_GpGp_exponential_anisotropic3D_alt', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn exponential_anisotropic3D_alt Derivatives of anisotropic Matern covariance
d_exponential_anisotropic3D_alt <- function(covparms, locs) {
    .Call('_GpGp_d_exponential_anisotropic3D_alt', PACKAGE = 'GpGp', covparms, locs)
}

#' Matern covariance function, different range parameter for each dimension
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range_1, ..., range_d, smoothness, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs is a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range_1, ..., range_d, smoothness, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range_1, ..., range_d, smoothness, nugget).
#' The covariance function is parameterized as
#' \deqn{ M(x,y) = \sigma^2 2^{1-\nu}/\Gamma(\nu) (|| D^{-1}(x - y) || )^\nu K_\nu(|| D^{-1}(x - y) || ) }
#' where D is a diagonal matrix with (range_1, ..., range_d) on the diagonals.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern_scaledim <- function(covparms, locs) {
    .Call('_GpGp_matern_scaledim', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern_scaledim Derivatives with respect to parameters
d_matern_scaledim <- function(covparms, locs) {
    .Call('_GpGp_d_matern_scaledim', PACKAGE = 'GpGp', covparms, locs)
}

#' Exponential covariance function, different range parameter for each dimension
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range_1, ..., range_d, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs is a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range_1, ..., range_d, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range_1, ..., range_d, nugget).
#' The covariance function is parameterized as
#' \deqn{ M(x,y) = \sigma^2 exp( - || D^{-1}(x - y) || ) }
#' where D is a diagonal matrix with (range_1, ..., range_d) on the diagonals.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
exponential_scaledim <- function(covparms, locs) {
    .Call('_GpGp_exponential_scaledim', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn exponential_scaledim Derivatives with respect to parameters
d_exponential_scaledim <- function(covparms, locs) {
    .Call('_GpGp_d_exponential_scaledim', PACKAGE = 'GpGp', covparms, locs)
}

#' Spatial-Temporal Matern covariance function
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range_1, range_2, smoothness, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d+1} columns.
#' Each row of locs is a point in R^{d+1}. The first \code{d} columns
#' should contain the spatial coordinates. The last column contains the times.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range_1, range_2, smoothness, nugget). range_1 is the
#' spatial range, and range_2 is the temporal range.
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range_1, range_2, smoothness, nugget).
#' The covariance function is parameterized as
#' \deqn{ M(x,y) = \sigma^2 2^{1-\nu}/\Gamma(\nu) (|| D^{-1}(x - y) || )^\nu K_\nu(|| D^{-1}(x - y) || ) }
#' where D is a diagonal matrix with (range_1, ..., range_1, range_2) on the diagonals.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern_spacetime <- function(covparms, locs) {
    .Call('_GpGp_matern_spacetime', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern_spacetime Derivatives with respect to parameters
d_matern_spacetime <- function(covparms, locs) {
    .Call('_GpGp_d_matern_spacetime', PACKAGE = 'GpGp', covparms, locs)
}

#' Spatial-Temporal exponential covariance function
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range_1, range_2, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d+1} columns.
#' Each row of locs is a point in R^{d+1}. The first \code{d} columns
#' should contain the spatial coordinates. The last column contains the times.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range_1, range_2, nugget). range_1 is the
#' spatial range, and range_2 is the temporal range.
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range_1, range_2, nugget).
#' The covariance function is parameterized as
#' \deqn{ M(x,y) = \sigma^2 exp( - || D^{-1}(x - y) || ) }
#' where D is a diagonal matrix with (range_1, ..., range_1, range_2) on the diagonals.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
exponential_spacetime <- function(covparms, locs) {
    .Call('_GpGp_exponential_spacetime', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn exponential_spacetime Derivatives with respect to parameters
d_exponential_spacetime <- function(covparms, locs) {
    .Call('_GpGp_d_exponential_spacetime', PACKAGE = 'GpGp', covparms, locs)
}

#' Isotropic Matern covariance function on sphere
#'
#' From a matrix of longitudes and latitudes and a vector covariance parameters of the form
#' (variance, range, smoothness, nugget), return the square matrix of
#' all pairwise covariances.
#' @param lonlat A matrix with \code{n} rows and one column with longitudes in (-180,180)
#' and one column of latitudes in (-90,90).
#' Each row of lonlat describes a point on the sphere.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range, smoothness, nugget). Range parameter assumes that
#' the sphere has radius 1 (units are radians).
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{lonlat[i,]} and
#' \code{lonlat[j,]}.
#' @section Matern on Sphere Domain:
#' The function first calculates the (x,y,z) 3D coordinates, and then inputs
#' the resulting locations into \code{matern_isotropic}. This means that we construct
#' covariances on the sphere by embedding the sphere in a 3D space. There has been some
#' concern expressed in the literature that such embeddings may produce distortions.
#' The source and nature of such distortions has never been articulated,
#' and to date, no such distortions have been documented. Guinness and
#' Fuentes (2016) argue that 3D embeddings produce reasonable models for data on spheres.
matern_sphere <- function(covparms, lonlat) {
    .Call('_GpGp_matern_sphere', PACKAGE = 'GpGp', covparms, lonlat)
}

#' @describeIn matern_sphere Derivatives with respect to parameters
d_matern_sphere <- function(covparms, lonlat) {
    .Call('_GpGp_d_matern_sphere', PACKAGE = 'GpGp', covparms, lonlat)
}

#' Isotropic exponential covariance function on sphere
#'
#' From a matrix of longitudes and latitudes and a vector covariance parameters of the form
#' (variance, range, nugget), return the square matrix of
#' all pairwise covariances.
#' @param lonlat A matrix with \code{n} rows and one column with longitudes in (-180,180)
#' and one column of latitudes in (-90,90).
#' Each row of lonlat describes a point on the sphere.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range, nugget). Range parameter assumes that
#' the sphere has radius 1 (units are radians).
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{lonlat[i,]} and
#' \code{lonlat[j,]}.
#' @section Covariances on spheres:
#' The function first calculates the (x,y,z) 3D coordinates, and then inputs
#' the resulting locations into \code{exponential_isotropic}. This means that we construct
#' covariances on the sphere by embedding the sphere in a 3D space. There has been some
#' concern expressed in the literature that such embeddings may produce distortions.
#' The source and nature of such distortions has never been articulated,
#' and to date, no such distortions have been documented. Guinness and
#' Fuentes (2016) argue that 3D embeddings produce reasonable models for data on spheres.
exponential_sphere <- function(covparms, lonlat) {
    .Call('_GpGp_exponential_sphere', PACKAGE = 'GpGp', covparms, lonlat)
}

#' @describeIn exponential_sphere Derivatives with respect to parameters
d_exponential_sphere <- function(covparms, lonlat) {
    .Call('_GpGp_d_exponential_sphere', PACKAGE = 'GpGp', covparms, lonlat)
}

#' Matern covariance function on sphere x time
#'
#' From a matrix of longitudes, latitudes, and times, and a vector covariance parameters of the form
#' (variance, range_1, range_2, smoothness, nugget), return the square matrix of
#' all pairwise covariances.
#' @param lonlattime A matrix with \code{n} rows and three columns: longitudes in (-180,180),
#' latitudes in (-90,90), and times.
#' Each row of lonlattime describes a point on the sphere x time.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range_1, range_2, smoothness, nugget), where range_1 is a 
#' spatial range (assuming sphere of radius 1), and range_2 is a temporal range.
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{lonlattime[i,]} and
#' \code{lonlattime[j,]}.
#' @section Covariances on spheres:
#' The function first calculates the (x,y,z) 3D coordinates, and then inputs
#' the resulting locations into \code{matern_spacetime}. This means that we construct
#' covariances on the sphere by embedding the sphere in a 3D space. There has been some
#' concern expressed in the literature that such embeddings may produce distortions.
#' The source and nature of such distortions has never been articulated,
#' and to date, no such distortions have been documented. Guinness and
#' Fuentes (2016) argue that 3D embeddings produce reasonable models for data on spheres.
matern_spheretime <- function(covparms, lonlattime) {
    .Call('_GpGp_matern_spheretime', PACKAGE = 'GpGp', covparms, lonlattime)
}

#' @describeIn matern_spheretime Derivatives with respect to parameters
d_matern_spheretime <- function(covparms, lonlattime) {
    .Call('_GpGp_d_matern_spheretime', PACKAGE = 'GpGp', covparms, lonlattime)
}

#' Exponential covariance function on sphere x time
#'
#' From a matrix of longitudes, latitudes, and times, and a vector covariance parameters of the form
#' (variance, range_1, range_2, nugget), return the square matrix of
#' all pairwise covariances.
#' @param lonlattime A matrix with \code{n} rows and three columns: longitudes in (-180,180),
#' latitudes in (-90,90), and times.
#' Each row of lonlattime describes a point on the sphere x time.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range_1, range_2, nugget), where range_1 is a 
#' spatial range (assuming sphere of radius 1), and range_2 is a temporal range.
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{lonlattime[i,]} and
#' \code{lonlattime[j,]}.
#' @section Covariances on spheres:
#' The function first calculates the (x,y,z) 3D coordinates, and then inputs
#' the resulting locations into \code{exponential_spacetime}. This means that we construct
#' covariances on the sphere by embedding the sphere in a 3D space. There has been some
#' concern expressed in the literature that such embeddings may produce distortions.
#' The source and nature of such distortions has never been articulated,
#' and to date, no such distortions have been documented. Guinness and
#' Fuentes (2016) argue that 3D embeddings produce reasonable models for data on spheres.
exponential_spheretime <- function(covparms, lonlattime) {
    .Call('_GpGp_exponential_spheretime', PACKAGE = 'GpGp', covparms, lonlattime)
}

#' @describeIn exponential_spheretime Derivatives with respect to parameters.
d_exponential_spheretime <- function(covparms, lonlattime) {
    .Call('_GpGp_d_exponential_spheretime', PACKAGE = 'GpGp', covparms, lonlattime)
}

#' Deformed Matern covariance function on sphere
#'
#' From a matrix of longitudes and latitudes and a vector covariance parameters of the form
#' (variance, range, smoothness, nugget, <5 warping parameters>), return the square matrix of
#' all pairwise covariances.
#' @param lonlat A matrix with \code{n} rows and one column with longitudes in (-180,180)
#' and one column of latitudes in (-90,90).
#' Each row of lonlat describes a point on the sphere.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range, smoothness, nugget, <5 warping parameters>). 
#' Range parameter assumes that the sphere has radius 1 (units are radians).
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{lonlat[i,]} and
#' \code{lonlat[j,]}.
#' @section Warpings:
#' The function first calculates the (x,y,z) 3D coordinates, and then "warps"
#' the locations to \eqn{(x,y,z) + \Phi(x,y,z)}, where \eqn{\Phi} is a warping
#' function composed of gradients of spherical harmonic functions of degree 2.
#' See Guinness (2019, "Gaussian Process Learning via Fisher Scoring of 
#' Vecchia's Approximation") for details.
#' The warped locations are input into \code{matern_isotropic}. 
matern_sphere_warp <- function(covparms, lonlat) {
    .Call('_GpGp_matern_sphere_warp', PACKAGE = 'GpGp', covparms, lonlat)
}

#' @describeIn matern_sphere_warp Derivatives with respect to parameters.
d_matern_sphere_warp <- function(covparms, lonlat) {
    .Call('_GpGp_d_matern_sphere_warp', PACKAGE = 'GpGp', covparms, lonlat)
}

#' Deformed exponential covariance function on sphere
#'
#' From a matrix of longitudes and latitudes and a vector covariance parameters of the form
#' (variance, range, nugget, <5 warping parameters>), return the square matrix of
#' all pairwise covariances.
#' @param lonlat A matrix with \code{n} rows and one column with longitudes in (-180,180)
#' and one column of latitudes in (-90,90).
#' Each row of lonlat describes a point on the sphere.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range, nugget, <5 warping parameters>). 
#' Range parameter assumes that the sphere has radius 1 (units are radians).
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{lonlat[i,]} and
#' \code{lonlat[j,]}.
#' @section Warpings:
#' The function first calculates the (x,y,z) 3D coordinates, and then "warps"
#' the locations to \eqn{(x,y,z) + \Phi(x,y,z)}, where \eqn{\Phi} is a warping
#' function composed of gradients of spherical harmonic functions of degree 2.
#' See Guinness (2019, "Gaussian Process Learning via Fisher Scoring of 
#' Vecchia's Approximation") for details.
#' The warped locations are input into \code{exponential_isotropic}. 
exponential_sphere_warp <- function(covparms, lonlat) {
    .Call('_GpGp_exponential_sphere_warp', PACKAGE = 'GpGp', covparms, lonlat)
}

#' @describeIn exponential_sphere_warp Derivatives with respect to parameters
d_exponential_sphere_warp <- function(covparms, lonlat) {
    .Call('_GpGp_d_exponential_sphere_warp', PACKAGE = 'GpGp', covparms, lonlat)
}

#' Deformed Matern covariance function on sphere
#'
#' From a matrix of longitudes, latitudes, times, and a vector covariance parameters of the form
#' (variance, range_1, range_2, smoothness, nugget, <5 warping parameters>), return the square matrix of
#' all pairwise covariances.
#' @param lonlattime A matrix with \code{n} rows and three columns: longitudes in (-180,180),
#' latitudes in (-90,90), and times.
#' Each row of lonlattime describes a point on the sphere x time.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range_1, range_2, smoothness, nugget, <5 warping parameters>). 
#' range_1 is a spatial range parameter that assumes that the sphere 
#' has radius 1 (units are radians). range_2 is a temporal range parameter.
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{lonlat[i,]} and
#' \code{lonlat[j,]}.
#' @section Warpings:
#' The function first calculates the (x,y,z) 3D coordinates, and then "warps"
#' the locations to \eqn{(x,y,z) + \Phi(x,y,z)}, where \eqn{\Phi} is a warping
#' function composed of gradients of spherical harmonic functions of degree 2.
#' See Guinness (2019, "Gaussian Process Learning via Fisher Scoring of 
#' Vecchia's Approximation") for details.
#' The warped locations are input into \code{matern_spacetime}. The function
#' does not do temporal warping.
matern_spheretime_warp <- function(covparms, lonlattime) {
    .Call('_GpGp_matern_spheretime_warp', PACKAGE = 'GpGp', covparms, lonlattime)
}

#' @describeIn matern_spheretime_warp Derivatives with respect to parameters
d_matern_spheretime_warp <- function(covparms, lonlattime) {
    .Call('_GpGp_d_matern_spheretime_warp', PACKAGE = 'GpGp', covparms, lonlattime)
}

#' Deformed exponential covariance function on sphere
#'
#' From a matrix of longitudes, latitudes, times, and a vector covariance parameters of the form
#' (variance, range_1, range_2, nugget, <5 warping parameters>), return the square matrix of
#' all pairwise covariances.
#' @param lonlattime A matrix with \code{n} rows and three columns: longitudes in (-180,180),
#' latitudes in (-90,90), and times.
#' Each row of lonlattime describes a point on the sphere x time.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range_1, range_2, nugget, <5 warping parameters>). 
#' range_1 is a spatial range parameter that assumes that the sphere 
#' has radius 1 (units are radians). range_2 is a temporal range parameter.
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{lonlat[i,]} and
#' \code{lonlat[j,]}.
#' @section Warpings:
#' The function first calculates the (x,y,z) 3D coordinates, and then "warps"
#' the locations to \eqn{(x,y,z) + \Phi(x,y,z)}, where \eqn{\Phi} is a warping
#' function composed of gradients of spherical harmonic functions of degree 2.
#' See Guinness (2019, "Gaussian Process Learning via Fisher Scoring of 
#' Vecchia's Approximation") for details.
#' The warped locations are input into \code{exponential_spacetime}. The function
#' does not do temporal warping.
exponential_spheretime_warp <- function(covparms, lonlattime) {
    .Call('_GpGp_exponential_spheretime_warp', PACKAGE = 'GpGp', covparms, lonlattime)
}

#' @describeIn exponential_spheretime_warp Derivatives with respect to parameters
d_exponential_spheretime_warp <- function(covparms, lonlattime) {
    .Call('_GpGp_d_exponential_spheretime_warp', PACKAGE = 'GpGp', covparms, lonlattime)
}

#' Isotropic Matern covariance function, nonstationary variances
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range, smoothness, nugget, <nonstat variance parameters>), 
#' return the square matrix of all pairwise covariances.
#' @param Z A matrix with \code{n} rows and \code{2} columns for spatial
#' locations + \code{p} columns describing spatial basis functions.
#' Each row of locs gives a point in R^2 (two dimensions only!) + the value
#' of \code{p} spatial basis functions.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range, smoothness, nugget, <nonstat variance parameters>).
#' The number of nonstationary variance parameters should equal \code{p}.
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' This covariance function multiplies the isotropic Matern covariance
#' by a nonstationary variance function. The form of the covariance is
#' \deqn{ C(x,y) = exp( \phi(x) + \phi(y) ) M(x,y) }
#' where M(x,y) is the isotropic Matern covariance, and 
#' \deqn{ \phi(x) = c_1 \phi_1(x) + ... + c_p \phi_p(x) }
#' where \eqn{\phi_1,...,\phi_p} are the spatial basis functions
#' contained in the last \code{p} columns of \code{Z}, and 
#' \eqn{c_1,...,c_p} are the nonstationary variance parameters.
matern_nonstat_var <- function(covparms, Z) {
    .Call('_GpGp_matern_nonstat_var', PACKAGE = 'GpGp', covparms, Z)
}

#' @describeIn matern_nonstat_var Derivatives with respect to parameters
d_matern_nonstat_var <- function(covparms, Z) {
    .Call('_GpGp_d_matern_nonstat_var', PACKAGE = 'GpGp', covparms, Z)
}

#' Isotropic exponential covariance function, nonstationary variances
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range, nugget, <nonstat variance parameters>), 
#' return the square matrix of all pairwise covariances.
#' @param Z A matrix with \code{n} rows and \code{2} columns for spatial
#' locations + \code{p} columns describing spatial basis functions.
#' Each row of locs gives a point in R^2 (two dimensions only!) + the value
#' of \code{p} spatial basis functions.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range, nugget, <nonstat variance parameters>).
#' The number of nonstationary variance parameters should equal \code{p}.
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' This covariance function multiplies the isotropic exponential covariance
#' by a nonstationary variance function. The form of the covariance is
#' \deqn{ C(x,y) = exp( \phi(x) + \phi(y) ) M(x,y) }
#' where M(x,y) is the isotropic exponential covariance, and 
#' \deqn{ \phi(x) = c_1 \phi_1(x) + ... + c_p \phi_p(x) }
#' where \eqn{\phi_1,...,\phi_p} are the spatial basis functions
#' contained in the last \code{p} columns of \code{Z}, and 
#' \eqn{c_1,...,c_p} are the nonstationary variance parameters.
exponential_nonstat_var <- function(covparms, Z) {
    .Call('_GpGp_exponential_nonstat_var', PACKAGE = 'GpGp', covparms, Z)
}

#' @describeIn exponential_nonstat_var Derivatives with respect to parameters
d_exponential_nonstat_var <- function(covparms, Z) {
    .Call('_GpGp_d_exponential_nonstat_var', PACKAGE = 'GpGp', covparms, Z)
}

#' Matern covariance function, smoothess = 1.5, different range parameter for each dimension
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range_1, ..., range_d, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs is a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range_1, ..., range_d, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range_1, ..., range_d, nugget).
#' The covariance function is parameterized as
#' \deqn{ M(x,y) = \sigma^2 (1 + || D^{-1}(x - y) || ) exp( - || D^{-1}(x - y) || ) }
#' where D is a diagonal matrix with (range_1, ..., range_d) on the diagonals.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern15_scaledim <- function(covparms, locs) {
    .Call('_GpGp_matern15_scaledim', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern15_scaledim Derivatives with respect to parameters
d_matern15_scaledim <- function(covparms, locs) {
    .Call('_GpGp_d_matern15_scaledim', PACKAGE = 'GpGp', covparms, locs)
}

#' Isotropic Matern covariance function, smoothness = 2.5
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs is a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range, nugget)
#' = \eqn{(\sigma^2,\alpha,\tau^2)}, and the covariance function is parameterized
#' as
#' \deqn{ M(x,y) = \sigma^2 (1 + || x - y ||/ \alpha + || x - y ||^2/3\alpha^2 ) exp( - || x - y ||/ \alpha )}
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern25_isotropic <- function(covparms, locs) {
    .Call('_GpGp_matern25_isotropic', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn exponential_isotropic Derivatives of isotropic
#' matern covariance function with smoothness 2.5
d_matern25_isotropic <- function(covparms, locs) {
    .Call('_GpGp_d_matern25_isotropic', PACKAGE = 'GpGp', covparms, locs)
}

#' Matern covariance function, smoothess = 2.5, different range parameter for each dimension
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range_1, ..., range_d, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs is a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range_1, ..., range_d, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range_1, ..., range_d, nugget).
#' The covariance function is parameterized as
#' \deqn{ M(x,y) = \sigma^2 (1 + || D^{-1}(x - y) || + || D^{-1}(x - y) ||^2/3.0) exp( - || D^{-1}(x - y) || ) }
#' where D is a diagonal matrix with (range_1, ..., range_d) on the diagonals.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern25_scaledim <- function(covparms, locs) {
    .Call('_GpGp_matern25_scaledim', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern25_scaledim Derivatives with respect to parameters
d_matern25_scaledim <- function(covparms, locs) {
    .Call('_GpGp_d_matern25_scaledim', PACKAGE = 'GpGp', covparms, locs)
}

#' Isotropic Matern covariance function, smoothness = 3.5
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs is a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range, nugget)
#' = \eqn{(\sigma^2,\alpha,\tau^2)}, and the covariance function is parameterized
#' as
#' \deqn{ M(x,y) = \sigma^2 ( \sum_{j=0}^3 c_j || x - y ||^j/ \alpha^j ) exp( - || x - y ||/ \alpha )}
#' where c_0 = 1, c_1 = 1, c_2 = 2/5, c_3 = 1/15. 
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern35_isotropic <- function(covparms, locs) {
    .Call('_GpGp_matern35_isotropic', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern35_isotropic Derivatives of isotropic
#' matern covariance function with smoothness 3.5
d_matern35_isotropic <- function(covparms, locs) {
    .Call('_GpGp_d_matern35_isotropic', PACKAGE = 'GpGp', covparms, locs)
}

#' Matern covariance function, smoothess = 3.5, different range parameter for each dimension
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range_1, ..., range_d, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs is a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range_1, ..., range_d, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range_1, ..., range_d, nugget).
#' The covariance function is parameterized as
#' \deqn{ M(x,y) = \sigma^2 ( \sum_{j=0}^3 c_j || D^{-1}(x - y) ||^j  ) exp( - || D^{-1}(x - y) || ) }
#' where c_0 = 1, c_1 = 1, c_2 = 2/5, c_3 = 1/15.
#' where D is a diagonal matrix with (range_1, ..., range_d) on the diagonals.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern35_scaledim <- function(covparms, locs) {
    .Call('_GpGp_matern35_scaledim', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern35_scaledim Derivatives with respect to parameters
d_matern35_scaledim <- function(covparms, locs) {
    .Call('_GpGp_d_matern35_scaledim', PACKAGE = 'GpGp', covparms, locs)
}

#' Isotropic Matern covariance function, smoothness = 4.5
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs is a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range, nugget)
#' = \eqn{(\sigma^2,\alpha,\tau^2)}, and the covariance function is parameterized
#' as
#' \deqn{ M(x,y) = \sigma^2 ( \sum_{j=0}^4 c_j || x - y ||^j/ \alpha^j ) exp( - || x - y ||/ \alpha )}
#' where c_0 = 1, c_1 = 1, c_2 = 3/7, c_3 = 2/21, c_4 = 1/105. 
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern45_isotropic <- function(covparms, locs) {
    .Call('_GpGp_matern45_isotropic', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern35_isotropic Derivatives of isotropic
#' matern covariance function with smoothness 3.5
d_matern45_isotropic <- function(covparms, locs) {
    .Call('_GpGp_d_matern45_isotropic', PACKAGE = 'GpGp', covparms, locs)
}

#' Matern covariance function, smoothess = 3.5, different range parameter for each dimension
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range_1, ..., range_d, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs is a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range_1, ..., range_d, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range_1, ..., range_d, nugget).
#' The covariance function is parameterized as
#' \deqn{ M(x,y) = \sigma^2 ( \sum_{j=0}^4 c_j || D^{-1}(x - y) ||^j  ) exp( - || D^{-1}(x - y) || ) }
#' where c_0 = 1, c_1 = 1, c_2 = 3/7, c_3 = 2/21, c_4 = 1/105. 
#' where D is a diagonal matrix with (range_1, ..., range_d) on the diagonals.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern45_scaledim <- function(covparms, locs) {
    .Call('_GpGp_matern45_scaledim', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern35_scaledim Derivatives with respect to parameters
d_matern45_scaledim <- function(covparms, locs) {
    .Call('_GpGp_d_matern45_scaledim', PACKAGE = 'GpGp', covparms, locs)
}

#' Isotropic Matern covariance function with random effects for categories
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, range, smoothness, category variance, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs gives a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, range, smoothness, category variance, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range, smoothness, category variance, nugget)
#' = \eqn{(\sigma^2,\alpha,\nu,c^2,\tau^2)}, and the covariance function is parameterized
#' as
#' \deqn{ M(x,y) = \sigma^2 2^{1-\nu}/\Gamma(\nu) (|| x - y ||/\alpha )^\nu K_\nu(|| x - y ||/\alpha ) }
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' The category variance \eqn{c^2} is added if two observation from same category
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern_categorical <- function(covparms, locs) {
    .Call('_GpGp_matern_categorical', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern_categorical Derivatives of isotropic Matern covariance
d_matern_categorical <- function(covparms, locs) {
    .Call('_GpGp_d_matern_categorical', PACKAGE = 'GpGp', covparms, locs)
}

#' Space-Time Matern covariance function with random effects for categories
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, spatial range, temporal range, smoothness, category, nugget), return the square matrix of
#' all pairwise covariances.
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs gives a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, spatial range, temporal range, smoothness, category, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range, smoothness, category, nugget)
#' = \eqn{(\sigma^2,\alpha_1,\alpha_2,\nu,c^2,\tau^2)}, and the covariance function is parameterized
#' as
#' \deqn{ d = ( || x - y ||^2/\alpha_1 + |s-t|^2/\alpha_2^2 )^{1/2} }
#' \deqn{ M(x,y) = \sigma^2 2^{1-\nu}/\Gamma(\nu) (d)^\nu K_\nu(d) }
#' (x,s) and (y,t) are the space-time locations of a pair of observations.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' The category variance \eqn{c^2} is added if two observation from same category
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern_spacetime_categorical <- function(covparms, locs) {
    .Call('_GpGp_matern_spacetime_categorical', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern_spacetime_categorical Derivatives of isotropic Matern covariance
d_matern_spacetime_categorical <- function(covparms, locs) {
    .Call('_GpGp_d_matern_spacetime_categorical', PACKAGE = 'GpGp', covparms, locs)
}

#' Space-Time Matern covariance function with local random effects for categories
#'
#' From a matrix of locations and covariance parameters of the form
#' (variance, spatial range, temporal range, smoothness, cat variance, cat spatial range, cat temporal range, cat smoothness, nugget),
#' return the square matrix of
#' all pairwise covariances.
#' This is the covariance for the following model for data from cateogory k
#' \deqn{ Y_k(x_i,t_i) = Z_0(x_i,t_i) + Z_k(x_i,t_i) + e_i }
#' where Z_0 is Matern with parameters (variance,spatial range,temporal range,smoothness)
#' and Z_1,...,Z_K are independent Materns with parameters
#' (cat variance, cat spatial range, cat temporal range, cat smoothness),
#' and e_1, ..., e_n are independent normals with variance (variance * nugget)
#' @param locs A matrix with \code{n} rows and \code{d} columns.
#' Each row of locs gives a point in R^d.
#' @param covparms A vector with covariance parameters
#' in the form (variance, spatial range, temporal range, smoothness, category, nugget)
#' @return A matrix with \code{n} rows and \code{n} columns, with the i,j entry
#' containing the covariance between observations at \code{locs[i,]} and
#' \code{locs[j,]}.
#' @section Parameterization:
#' The covariance parameter vector is (variance, range, smoothness, category, nugget)
#' = \eqn{(\sigma^2,\alpha_1,\alpha_2,\nu,c^2,\tau^2)}, and the covariance function is parameterized
#' as
#' \deqn{ d = ( || x - y ||^2/\alpha_1 + |s-t|^2/\alpha_2^2 )^{1/2} }
#' \deqn{ M(x,y) = \sigma^2 2^{1-\nu}/\Gamma(\nu) (d)^\nu K_\nu(d) }
#' (x,s) and (y,t) are the space-time locations of a pair of observations.
#' The nugget value \eqn{ \sigma^2 \tau^2 } is added to the diagonal of the covariance matrix.
#' The category variance \eqn{c^2} is added if two observation from same category
#' NOTE: the nugget is \eqn{ \sigma^2 \tau^2 }, not \eqn{ \tau^2 }. 
matern_spacetime_categorical_local <- function(covparms, locs) {
    .Call('_GpGp_matern_spacetime_categorical_local', PACKAGE = 'GpGp', covparms, locs)
}

#' @describeIn matern_spacetime_categorical_local Derivatives of isotropic Matern covariance
d_matern_spacetime_categorical_local <- function(covparms, locs) {
    .Call('_GpGp_d_matern_spacetime_categorical_local', PACKAGE = 'GpGp', covparms, locs)
}

#' Multiply approximate inverse Cholesky by a vector
#'
#' Vecchia's approximation implies a sparse approximation to the
#' inverse Cholesky factor of the covariance matrix. This function
#' returns the result of multiplying that matrix by a vector.
#' @param Linv Entries of the sparse inverse Cholesky factor,
#' usually the output from \code{\link{vecchia_Linv}}.
#' @param z the vector to be multiplied
#' @inheritParams vecchia_meanzero_loglik
#' @return the product of the sparse inverse Cholesky factor with a vector
#' @examples
#' n <- 2000
#' locs <- matrix( runif(2*n), n, 2 )
#' covparms <- c(2, 0.2, 0.75, 0.1)
#' ord <- order_maxmin(locs)
#' NNarray <- find_ordered_nn(locs,20)
#' Linv <- vecchia_Linv( covparms, "matern_isotropic", locs, NNarray )
#' z1 <- rnorm(n)
#' y <- fast_Gp_sim_Linv(Linv,NNarray,z1)
#' z2 <- Linv_mult(Linv, y, NNarray)
#' print( sum( (z1-z2)^2 ) )
#' @export
Linv_mult <- function(Linv, z, NNarray) {
    .Call('_GpGp_Linv_mult', PACKAGE = 'GpGp', Linv, z, NNarray)
}

#' Multiply approximate Cholesky by a vector
#'
#' Vecchia's approximation implies a sparse approximation to the
#' inverse Cholesky factor of the covariance matrix. This function
#' returns the result of multiplying the inverse of that matrix by a vector
#' (i.e. an approximation to the Cholesky factor).
#' @param Linv Entries of the sparse inverse Cholesky factor,
#' usually the output from \code{\link{vecchia_Linv}}.
#' @param z the vector to be multiplied
#' @inheritParams vecchia_meanzero_loglik
#' @return the product of the Cholesky factor with a vector
#' @examples
#' n <- 2000
#' locs <- matrix( runif(2*n), n, 2 )
#' covparms <- c(2, 0.2, 0.75, 0.1)
#' ord <- order_maxmin(locs)
#' NNarray <- find_ordered_nn(locs,20)
#' Linv <- vecchia_Linv( covparms, "matern_isotropic", locs, NNarray )
#' z <- rnorm(n)
#' y1 <- fast_Gp_sim_Linv(Linv,NNarray,z)
#' y2 <- L_mult(Linv, z, NNarray)
#' print( sum( (y1-y2)^2 ) )
#' @export
L_mult <- function(Linv, z, NNarray) {
    .Call('_GpGp_L_mult', PACKAGE = 'GpGp', Linv, z, NNarray)
}

#' Multiply transpose of approximate inverse Cholesky by a vector
#'
#' Vecchia's approximation implies a sparse approximation to the
#' inverse Cholesky factor of the covariance matrix. This function
#' returns the result of multiplying the transpose of that matrix by a vector.
#' @param Linv Entries of the sparse inverse Cholesky factor,
#' usually the output from \code{\link{vecchia_Linv}}.
#' @param z the vector to be multiplied
#' @inheritParams vecchia_meanzero_loglik
#' @return the product of the transpose of the 
#' sparse inverse Cholesky factor with a vector
#' @examples
#' n <- 2000
#' locs <- matrix( runif(2*n), n, 2 )
#' covparms <- c(2, 0.2, 0.75, 0.1)
#' NNarray <- find_ordered_nn(locs,20)
#' Linv <- vecchia_Linv( covparms, "matern_isotropic", locs, NNarray )
#' z1 <- rnorm(n)
#' z2 <- Linv_t_mult(Linv, z1, NNarray)
#' @export
Linv_t_mult <- function(Linv, z, NNarray) {
    .Call('_GpGp_Linv_t_mult', PACKAGE = 'GpGp', Linv, z, NNarray)
}

#' Multiply transpose of approximate Cholesky by a vector
#'
#' Vecchia's approximation implies a sparse approximation to the
#' inverse Cholesky factor of the covariance matrix. This function
#' returns the result of multiplying the transpose of the
#' inverse of that matrix by a vector
#' (i.e. an approximation to the transpose of the Cholesky factor).
#' @param Linv Entries of the sparse inverse Cholesky factor,
#' usually the output from \code{\link{vecchia_Linv}}.
#' @param z the vector to be multiplied
#' @inheritParams vecchia_meanzero_loglik
#' @return the product of the transpose of the Cholesky factor with a vector
#' @examples
#' n <- 2000
#' locs <- matrix( runif(2*n), n, 2 )
#' covparms <- c(2, 0.2, 0.75, 0.1)
#' NNarray <- find_ordered_nn(locs,20)
#' Linv <- vecchia_Linv( covparms, "matern_isotropic", locs, NNarray )
#' z1 <- rnorm(n)
#' z2 <- L_t_mult(Linv, z1, NNarray)
#' @export
L_t_mult <- function(Linv, z, NNarray) {
    .Call('_GpGp_L_t_mult', PACKAGE = 'GpGp', Linv, z, NNarray)
}

#' Entries of inverse Cholesky approximation
#' 
#' This function returns the entries of the inverse Cholesky
#' factor of the covariance matrix implied by Vecchia's approximation.
#' For return matrix \code{Linv}, \code{Linv[i,]} contains 
#' the non-zero entries of row \code{i} of
#' the inverse Cholesky matrix. The columns of the non-zero entries
#' are specified in \code{NNarray[i,]}.
#' @inheritParams vecchia_meanzero_loglik
#' @param start_ind Compute entries of Linv only for rows \code{start_ind}
#' until the last row. 
#' @return matrix containing entries of inverse Cholesky
#' @examples
#' n1 <- 40
#' n2 <- 40
#' n <- n1*n2
#' locs <- as.matrix( expand.grid( (1:n1)/n1, (1:n2)/n2 ) )
#' covparms <- c(2, 0.2, 0.75, 0)
#' NNarray <- find_ordered_nn(locs,20)
#' Linv <- vecchia_Linv(covparms, "matern_isotropic", locs, NNarray)
#' @export
vecchia_Linv <- function(covparms, covfun_name, locs, NNarray, start_ind = 1L) {
    .Call('_GpGp_vecchia_Linv', PACKAGE = 'GpGp', covparms, covfun_name, locs, NNarray, start_ind)
}

#' Vecchia's loglikelihood, gradient, and Fisher information
#'
#' This function returns Vecchia's (1988) approximation to the Gaussian
#' loglikelihood, profiling out the regression coefficients, and returning
#' the gradient and Fisher information. 
#' Vecchia's approximation modifies the ordered conditional
#' specification of the joint density; rather than each term in the product
#' conditioning on all previous observations, each term conditions on
#' a small subset of previous observations.
#' @inheritParams vecchia_meanzero_loglik
#' @param X Design matrix of covariates. Row \code{i} of \code{X} contains
#' the covariates for the observation at row \code{i} of \code{locs}.
#' @return A list containing 
#' \itemize{
#'     \item \code{loglik}: the loglikelihood
#'     \item \code{grad}: gradient with respect to covariance parameters
#'     \item \code{info}: Fisher information for covariance parameters
#'     \item \code{betahat}: profile likelihood estimate of regression coefs
#'     \item \code{betainfo}: information matrix for \code{betahat}.
#' }
#' The covariance matrix for \code{$betahat} is the inverse of \code{$betainfo}.
#' @examples
#' n1 <- 20
#' n2 <- 20
#' n <- n1*n2
#' locs <- as.matrix( expand.grid( (1:n1)/n1, (1:n2)/n2 ) )
#' X <- cbind(rep(1,n),locs[,2])
#' covparms <- c(2, 0.2, 0.75, 0)
#' y <- X %*% c(1,2) + fast_Gp_sim(covparms, "matern_isotropic", locs, 50 )
#' ord <- order_maxmin(locs)
#' NNarray <- find_ordered_nn(locs,20)
#' #loglik <- vecchia_profbeta_loglik_grad_info( covparms, "matern_isotropic", 
#' #    y, X, locs, NNarray )
#' @export
vecchia_profbeta_loglik_grad_info <- function(covparms, covfun_name, y, X, locs, NNarray) {
    .Call('_GpGp_vecchia_profbeta_loglik_grad_info', PACKAGE = 'GpGp', covparms, covfun_name, y, X, locs, NNarray)
}

#' Vecchia's approximation to the Gaussian loglikelihood, with profiled 
#' regression coefficients.
#'
#' This function returns Vecchia's (1988) approximation to the Gaussian
#' loglikelihood, profiling out the regression coefficients. 
#' The approximation modifies the ordered conditional
#' specification of the joint density; rather than each term in the product
#' conditioning on all previous observations, each term conditions on
#' a small subset of previous observations.
#' @inheritParams vecchia_meanzero_loglik
#' @param X Design matrix of covariates. Row \code{i} of \code{X} contains
#' the covariates for the observation at row \code{i} of \code{locs}.
#' @return a list containing
#' \itemize{
#'  \item \code{loglik}: the loglikelihood
#'  \item \code{betahat}: profile likelihood estimate of regression coefficients
#'  \item \code{betainfo}: information matrix for \code{betahat}.
#' }
#' The covariance
#' matrix for \code{$betahat} is the inverse of \code{$betainfo}.
#' @examples
#' n1 <- 20
#' n2 <- 20
#' n <- n1*n2
#' locs <- as.matrix( expand.grid( (1:n1)/n1, (1:n2)/n2 ) )
#' X <- cbind(rep(1,n),locs[,2])
#' covparms <- c(2, 0.2, 0.75, 0)
#' y <- X %*% c(1,2) + fast_Gp_sim(covparms, "matern_isotropic", locs, 50 )
#' ord <- order_maxmin(locs)
#' NNarray <- find_ordered_nn(locs,20)
#' #loglik <- vecchia_profbeta_loglik( covparms, "matern_isotropic", y, X, locs, NNarray )
#' @export
vecchia_profbeta_loglik <- function(covparms, covfun_name, y, X, locs, NNarray) {
    .Call('_GpGp_vecchia_profbeta_loglik', PACKAGE = 'GpGp', covparms, covfun_name, y, X, locs, NNarray)
}

#' Vecchia's approximation to the Gaussian loglikelihood, zero mean
#'
#' This function returns Vecchia's (1988) approximation to the Gaussian
#' loglikelihood. The approximation modifies the ordered conditional
#' specification of the joint density; rather than each term in the product
#' conditioning on all previous observations, each term conditions on
#' a small subset of previous observations.
#' @param covparms A vector of covariance parameters appropriate
#' for the specified covariance function
#' @param covfun_name See \code{\link{GpGp}} for information about covariance
#' functions.
#' @param y vector of response values
#' @param locs matrix of locations. Row \code{i} of \code{locs} specifies the location
#' of element \code{i} of \code{y}, and so the length of \code{y} should equal
#' the number of rows of \code{locs}.
#' @param NNarray A matrix of indices, usually the output from \code{\link{find_ordered_nn}}.
#' Row \code{i} contains the indices
#' of the observations that observation \code{i} conditions on. By convention,
#' the first element of row \code{i} is \code{i}.
#' @return a list containing
#' \itemize{
#'  \item \code{loglik}: the loglikelihood
#' }
#' @examples
#' n1 <- 20
#' n2 <- 20
#' n <- n1*n2
#' locs <- as.matrix( expand.grid( (1:n1)/n1, (1:n2)/n2 ) )
#' covparms <- c(2, 0.2, 0.75, 0)
#' y <- fast_Gp_sim(covparms, "matern_isotropic", locs, 50 )
#' ord <- order_maxmin(locs)
#' NNarray <- find_ordered_nn(locs,20)
#' #loglik <- vecchia_meanzero_loglik( covparms, "matern_isotropic", y, locs, NNarray )
#' @export
vecchia_meanzero_loglik <- function(covparms, covfun_name, y, locs, NNarray) {
    .Call('_GpGp_vecchia_meanzero_loglik', PACKAGE = 'GpGp', covparms, covfun_name, y, locs, NNarray)
}

#' Grouped Vecchia loglikelihood, gradient, Fisher information
#'
#' This function returns a grouped version (Guinness, 2018) of Vecchia's (1988) 
#' approximation to the Gaussian
#' loglikelihood, the gradient, and Fisher information, 
#' and the profile likelihood estimate of the regression
#' coefficients. The approximation modifies the ordered conditional
#' specification of the joint density; rather than each term in the product
#' conditioning on all previous observations, each term conditions on
#' a small subset of previous observations.
#' @inheritParams vecchia_grouped_meanzero_loglik
#' @param X Design matrix of covariates. Row \code{i} of \code{X} contains
#' the covariates for the observation at row \code{i} of \code{locs}.
#' @return a list containing
#' \itemize{
#'     \item \code{loglik}: the loglikelihood
#'     \item \code{grad}: gradient with respect to covariance parameters
#'     \item \code{info}: Fisher information for covariance parameters
#'     \item \code{betahat}: profile likelihood estimate of regression coefs
#'     \item \code{betainfo}: information matrix for \code{betahat}.
#' }
#' The covariance
#' matrix for \code{$betahat} is the inverse of \code{$betainfo}.
#' @examples
#' n1 <- 20
#' n2 <- 20
#' n <- n1*n2
#' locs <- as.matrix( expand.grid( (1:n1)/n1, (1:n2)/n2 ) )
#' X <- cbind(rep(1,n),locs[,2])
#' covparms <- c(2, 0.2, 0.75, 0)
#' y <- fast_Gp_sim(covparms, "matern_isotropic", locs, 50 )
#' ord <- order_maxmin(locs)
#' NNarray <- find_ordered_nn(locs,20)
#' NNlist <- group_obs(NNarray)
#' #loglik <- vecchia_grouped_profbeta_loglik_grad_info( 
#' #    covparms, "matern_isotropic", y, X, locs, NNlist )
#' @export
vecchia_grouped_profbeta_loglik_grad_info <- function(covparms, covfun_name, y, X, locs, NNlist) {
    .Call('_GpGp_vecchia_grouped_profbeta_loglik_grad_info', PACKAGE = 'GpGp', covparms, covfun_name, y, X, locs, NNlist)
}

#' Grouped Vecchia approximation, profiled regression coefficients
#'
#' This function returns a grouped version (Guinness, 2018) of Vecchia's (1988) 
#' approximation to the Gaussian
#' loglikelihood and the profile likelihood estimate of the regression
#' coefficients. The approximation modifies the ordered conditional
#' specification of the joint density; rather than each term in the product
#' conditioning on all previous observations, each term conditions on
#' a small subset of previous observations.
#' @inheritParams vecchia_grouped_meanzero_loglik
#' @param X Design matrix of covariates. Row \code{i} of \code{X} contains
#' the covariates for the observation at row \code{i} of \code{locs}.
#' @return a list containing
#' \itemize{
#'  \item \code{loglik}: the loglikelihood
#'  \item \code{betahat}: profile likelihood estimate of regression coefficients
#'  \item \code{betainfo}: information matrix for \code{betahat}.
#' }
#' The covariance
#' matrix for \code{$betahat} is the inverse of \code{$betainfo}.
#' @examples
#' n1 <- 20
#' n2 <- 20
#' n <- n1*n2
#' locs <- as.matrix( expand.grid( (1:n1)/n1, (1:n2)/n2 ) )
#' X <- cbind(rep(1,n),locs[,2])
#' covparms <- c(2, 0.2, 0.75, 0)
#' y <- fast_Gp_sim(covparms, "matern_isotropic", locs, 50 )
#' ord <- order_maxmin(locs)
#' NNarray <- find_ordered_nn(locs,20)
#' NNlist <- group_obs(NNarray)
#' #loglik <- vecchia_grouped_profbeta_loglik( 
#' #    covparms, "matern_isotropic", y, X, locs, NNlist )
#' @export
vecchia_grouped_profbeta_loglik <- function(covparms, covfun_name, y, X, locs, NNlist) {
    .Call('_GpGp_vecchia_grouped_profbeta_loglik', PACKAGE = 'GpGp', covparms, covfun_name, y, X, locs, NNlist)
}

#' Grouped Vecchia approximation to the Gaussian loglikelihood, zero mean
#'
#' This function returns a grouped version (Guinness, 2018) of Vecchia's (1988) 
#' approximation to the Gaussian
#' loglikelihood. The approximation modifies the ordered conditional
#' specification of the joint density; rather than each term in the product
#' conditioning on all previous observations, each term conditions on
#' a small subset of previous observations.
#' @param covparms A vector of covariance parameters appropriate
#' for the specified covariance function
#' @param covfun_name See \code{\link{GpGp}} for information about covariance
#' functions.
#' @param y vector of response values
#' @param locs matrix of locations. Row \code{i} of \code{locs} specifies the location
#' of element \code{i} of \code{y}, and so the length of \code{y} should equal
#' the number of rows of \code{locs}.
#' @param NNlist A neighbor list object, the output from \code{\link{group_obs}}.
#' @return a list containing
#' \itemize{
#'  \item \code{loglik}: the loglikelihood
#' }
#' @examples
#' n1 <- 20
#' n2 <- 20
#' n <- n1*n2
#' locs <- as.matrix( expand.grid( (1:n1)/n1, (1:n2)/n2 ) )
#' covparms <- c(2, 0.2, 0.75, 0)
#' y <- fast_Gp_sim(covparms, "matern_isotropic", locs, 50 )
#' ord <- order_maxmin(locs)
#' NNarray <- find_ordered_nn(locs,20)
#' NNlist <- group_obs(NNarray)
#' #loglik <- vecchia_grouped_meanzero_loglik( covparms, "matern_isotropic", y, locs, NNlist )
#' @export
vecchia_grouped_meanzero_loglik <- function(covparms, covfun_name, y, locs, NNlist) {
    .Call('_GpGp_vecchia_grouped_meanzero_loglik', PACKAGE = 'GpGp', covparms, covfun_name, y, locs, NNlist)
}