R/RcppExports.R

In polykde: Polyspherical Kernel Density Estimation

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

#' @title Euler algorithms for polyspherical density ridge estimation
#'
#' @description Functions to perform density ridge estimation on the
#' polysphere \eqn{\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}}
#' through the Euler algorithm in standard, parallel, or block mode.
#'
#' @param x a matrix of size \code{c(nx, sum(d) + r)} with the starting points
#' for the Euler algorithm.
#' @inheritParams kde_polysph
#' @inheritParams proj_grad_kde_polysph
#' @param h_euler vector of size \code{r} with the advance steps in the Euler
#' method. Set internally as \code{h} if not provided.
#' @param N maximum number of Euler iterations. Defaults to \code{1e3}.
#' @param eps convergence tolerance. Defaults to \code{1e-5}.
#' @param keep_paths keep the Euler paths to the ridge? Defaults to
#' \code{FALSE}.
#' @param show_prog display a progress bar for \code{x}? Defaults to
#' \code{TRUE}.
#' @param show_prog_j display a progress bar for \code{N}? Defaults to
#' \code{FALSE}.
#' @param ind_blocks indexes of the blocks, a vector or length \code{r}.
#' @param ... further arguments passed to \code{\link{euler_ridge}}.
#' @param cores cores to use. Defaults to \code{1}.
#' @details \code{euler_ridge} is the main function to perform density ridge
#' estimation through the Euler algorithm from the starting values \code{x}
#' to initiate the ridge path. The function \code{euler_ridge_parallel}
#' parallelizes on the starting values \code{x}. The function
#' \code{euler_ridge_block} runs the Euler algorithm marginally in blocks
#' of spheres, instead of jointly in the whole polysphere. This function
#' requires that all the dimensions are the same.
#' @return The three functions return a list with the following fields:
#' \item{ridge_y}{a matrix of size \code{c(nx, sum(d) + r)} with the end
#' points of Euler algorithm defining the estimated ridge.}
#' \item{lamb_norm_y}{a matrix of size \code{c(nx, sum(d) + r)} with the
#' Hessian eigenvalues (largest to smallest) evaluated at end points.}
#' \item{log_dens_y}{a column vector of size \code{c(nx, 1)} with the
#' logarithm of the density at end points.}
#' \item{paths}{an array of size \code{c(nx, sum(d) + r, N + 1)} containing
#' the Euler paths.}
#' \item{start_x}{a matrix of size \code{c(nx, sum(d) + r)} with the starting
#' points for the Euler algorithm.}
#' \item{iter}{a column vector of size \code{c(nx, 1)} counting the iterations
#' required for each point.}
#' \item{conv}{a column vector of size \code{c(nx, 1)} with convergence flags.}
#' \item{d}{vector \code{d}.}
#' \item{h}{bandwidth used for the kernel density estimator.}
#' \item{error}{a column vector of size \code{c(nx, 1)} indicating if errors
#' were found for each path.}
#' @examples
#' ## Test on S^2 with a small circle trend
#'
#' # Sample
#' r <- 1
#' d <- 2
#' n <- 50
#' ind_dj <- comp_ind_dj(d = d)
#' set.seed(987204452)
#' X <- r_path_s2r(n = n, r = r, spiral = FALSE, Theta = cbind(c(1, 0, 0)),
#'                 sigma = 0.35)[, , 1]
#' col_X_alp <- viridis::viridis(n, alpha = 0.25)
#' col_X <- viridis::viridis(n)
#'
#' # Euler
#' h_rid <- 0.5
#' h_eu <- h_rid^2
#' N <- 30
#' eps <- 1e-6
#' Y <- euler_ridge(x = X, X = X, d = d, h = h_rid, h_euler = h_eu,
#'                  N = N, eps = eps, keep_paths = TRUE)
#' Y
#'
#' # Visualization
#' i <- N # Between 1 and N
#' sc3 <- scatterplot3d::scatterplot3d(Y$paths[, , 1], color = col_X_alp,
#'                                     pch = 19, xlim = c(-1, 1),
#'                                     ylim = c(-1, 1), zlim = c(-1, 1),
#'                                     xlab = "x", ylab = "y", zlab = "z")
#' sc3$points3d(rbind(Y$paths[, , i]), col = col_X, pch = 16, cex = 0.75)
#' for (k in seq_len(nrow(Y$paths))) {
#'
#'   sc3$points3d(t(Y$paths[k, , ]), col = col_X_alp[k], type = "l")
#'
#' }
#' @export
euler_ridge <- function(x, X, d, h, h_euler = as.numeric( c()), weights = as.numeric( c()), wrt_unif = FALSE, normalized = TRUE, norm_x = FALSE, norm_X = FALSE, kernel = 1L, kernel_type = 1L, k = 10.0, N = 1e3L, eps = 1e-5, keep_paths = FALSE, proj_alt = TRUE, fix_u1 = TRUE, sparse = FALSE, show_prog = TRUE, show_prog_j = FALSE) {
    .Call(`_polykde_euler_ridge`, x, X, d, h, h_euler, weights, wrt_unif, normalized, norm_x, norm_X, kernel, kernel_type, k, N, eps, keep_paths, proj_alt, fix_u1, sparse, show_prog, show_prog_j)
}

#' @title Gradient and Hessian of the polyspherical kernel density estimator
#'
#' @description Computes the gradient
#' \eqn{\mathsf{D}\hat{f}(\boldsymbol{x};\boldsymbol{h})} and Hessian matrix
#' \eqn{\mathsf{H}\hat{f}(\boldsymbol{x};\boldsymbol{h})} of the kernel density
#' estimator \eqn{\hat{f}(\boldsymbol{x};\boldsymbol{h})} on the
#' polysphere \eqn{\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}}.
#'
#' @inheritParams kde_polysph
#' @param projected compute the \emph{projected} gradient and Hessian that
#' accounts for the radial projection? Defaults to \code{TRUE}.
#' @param proj_alt alternative projection. Defaults to \code{TRUE}.
#' @param norm_grad_hess normalize the gradient and Hessian dividing by the
#' kernel density estimator? Defaults to \code{FALSE}.
#' @return A list with the following components:
#' \item{dens}{a column vector of size \code{c(nx, 1)} with the kernel
#' density estimator evaluated at \code{x}.}
#' \item{grad}{a matrix of size \code{c(nx, sum(d) + r)} with the gradient of
#' the kernel density estimator evaluated at \code{x}.}
#' \item{hess}{an array of size \code{c(nx, sum(d) + r, sum(d) + r)} with the
#' Hessian matrix of the kernel density estimator evaluated at \code{x}.}
#' @examples
#' # Simple check on (S^1)^2
#' n <- 3
#' d <- c(1, 1)
#' mu <- c(0, 1, 0, 1)
#' kappa <- c(5, 5)
#' h <- c(0.2, 0.2)
#' X <- r_vmf_polysph(n = n, d = d, mu = mu, kappa = kappa)
#' grh <- grad_hess_kde_polysph(x = X, X = X, d = d, h = h)
#' str(grh)
#' grh
#' @export
grad_hess_kde_polysph <- function(x, X, d, h, weights = as.numeric( c()), projected = TRUE, proj_alt = TRUE, norm_grad_hess = FALSE, log = FALSE, wrt_unif = FALSE, normalized = TRUE, norm_x = FALSE, norm_X = FALSE, kernel = 1L, kernel_type = 1L, k = 10.0) {
    .Call(`_polykde_grad_hess_kde_polysph`, x, X, d, h, weights, projected, proj_alt, norm_grad_hess, log, wrt_unif, normalized, norm_x, norm_X, kernel, kernel_type, k)
}

#' @title Projected gradient of the polyspherical kernel density estimator
#'
#' @description Computes the projected gradient
#' \eqn{\mathsf{D}_{(p-1)}\hat{f}(\boldsymbol{x};\boldsymbol{h})} of the
#' kernel density estimator \eqn{\hat{f}(\boldsymbol{x};\boldsymbol{h})} on the
#' polysphere \eqn{\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}},
#' where \eqn{p=\sum_{j=1}^r d_j+r} is the dimension of the ambient space.
#'
#' @inheritParams kde_polysph
#' @inheritParams grad_hess_kde_polysph
#' @param fix_u1 ensure the \eqn{u_1} vector is different from \eqn{x}?
#' Prevents the Euler algorithm to "surf the ridge". Defaults to \code{TRUE}.
#' @param sparse use a sparse eigendecomposition of the Hessian? Defaults to
#' \code{FALSE}.
#' @return A list with the following components:
#' \item{eta}{a matrix of size \code{c(nx, sum(d) + r)} with the
#' projected gradient evaluated at \code{x}.}
#' \item{u1}{a matrix of size \code{c(nx, sum(d) + r)} with the
#' first non-null Hessian eigenvector evaluated at \code{x}.}
#' \item{lamb_norm}{a matrix of size \code{c(nx, sum(d) + r)} with the
#' Hessian eigenvalues (largest to smallest) evaluated at \code{x}.}
#' @examples
#' # Simple check on (S^1)^2
#' n <- 3
#' d <- c(1, 1)
#' mu <- c(0, 1, 0, 1)
#' kappa <- c(5, 5)
#' h <- c(0.2, 0.2)
#' X <- r_vmf_polysph(n = n, d = d, mu = mu, kappa = kappa)
#' proj_grad_kde_polysph(x = X, X = X, d = d, h = h)
#' @export
proj_grad_kde_polysph <- function(x, X, d, h, weights = as.numeric( c()), wrt_unif = FALSE, normalized = TRUE, norm_x = FALSE, norm_X = FALSE, kernel = 1L, kernel_type = 1L, k = 10.0, proj_alt = TRUE, fix_u1 = TRUE, sparse = FALSE) {
    .Call(`_polykde_proj_grad_kde_polysph`, x, X, d, h, weights, wrt_unif, normalized, norm_x, norm_X, kernel, kernel_type, k, proj_alt, fix_u1, sparse)
}

#' @title Polyspherical kernel density estimator
#'
#' @description Computes the kernel density estimator for data on the
#' polysphere \eqn{\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}}.
#' Given a sample \eqn{\boldsymbol{X}_1,\ldots,\boldsymbol{X}_n}, this
#' estimator is
#' \deqn{\hat{f}(\boldsymbol{x};\boldsymbol{h})=\sum_{i=1}^n
#' L_{\boldsymbol{h}}(\boldsymbol{x},\boldsymbol{X}_i)}
#' for a kernel \eqn{L} and a vector of bandwidths \eqn{\boldsymbol{h}}.
#'
#' @param x a matrix of size \code{c(nx, sum(d) + r)} with the evaluation
#' points.
#' @param X a matrix of size \code{c(n, sum(d) + r)} with the sample.
#' @param d vector of size \code{r} with dimensions.
#' @param h vector of size \code{r} with bandwidths.
#' @param weights weights for each observation. If provided, a vector of size
#' \code{n} with the weights for multiplying each kernel. If not provided,
#' set internally to \code{rep(1 / n, n)}, which gives the standard estimator.
#' @param log compute the logarithm of the density? Defaults to \code{FALSE}.
#' @param wrt_unif flag to return a density with respect to the uniform
#' measure. If \code{FALSE} (default), the density is with respect to the
#' Lebesgue measure.
#' @param normalized flag to compute the normalizing constant of the kernel
#' and include it in the kernel density estimator. Defaults to \code{TRUE}.
#' @param intrinsic use the intrinsic distance, instead of the
#' extrinsic-chordal distance, in the kernel? Defaults to \code{FALSE}.
#' @param norm_x,norm_X ensure a normalization of the data? Defaults to
#' \code{FALSE}.
#' @param kernel kernel employed: \code{1} for von Mises--Fisher (default);
#' \code{2} for Epanechnikov; \code{3} for softplus.
#' @param kernel_type type of kernel employed: \code{1} for product kernel
#' (default); \code{2} for spherically symmetric kernel.
#' @param k softplus kernel parameter. Defaults to \code{10.0}.
#' @return A column vector of size \code{c(nx, 1)} with the evaluation of
#' kernel density estimator.
#' @examples
#' # Simple check on S^1 x S^2
#' n <- 1e3
#' d <- c(1, 2)
#' mu <- c(0, 1, 0, 0, 1)
#' kappa <- c(5, 5)
#' h <- c(0.2, 0.2)
#' X <- r_vmf_polysph(n = n, d = d, mu = mu, kappa = kappa)
#' kde_polysph(x = rbind(mu), X = X, d = d, h = h)
#' d_vmf_polysph(x = rbind(mu), d = d, mu = mu, kappa = kappa)
#' @export
kde_polysph <- function(x, X, d, h, weights = as.numeric( c()), log = FALSE, wrt_unif = FALSE, normalized = TRUE, intrinsic = FALSE, norm_x = FALSE, norm_X = FALSE, kernel = 1L, kernel_type = 1L, k = 10.0) {
    .Call(`_polykde_kde_polysph`, x, X, d, h, weights, log, wrt_unif, normalized, intrinsic, norm_x, norm_X, kernel, kernel_type, k)
}

#' @title Cross-validation for the polyspherical kernel density estimator
#'
#' @description Computes the logarithm of the cross-validated kernel density
#' estimator: \eqn{\log \hat{f}_{-i}(\boldsymbol{X}_i;\boldsymbol{h})},
#' \eqn{i = 1, \ldots, n.}
#'
#' @inheritParams kde_polysph
#' @param norm_X ensure a normalization of the data? Defaults to \code{FALSE}.
#' @return A column vector of size \code{c(n, 1)} with the evaluation of the
#' logarithm of the cross-validated kernel density estimator.
#' @examples
#' # Simple check on S^1 x S^2
#' n <- 5
#' d <- c(1, 2)
#' h <- c(0.2, 0.2)
#' X <- r_unif_polysph(n = n, d = d)
#' log_cv_kde_polysph(X = X, d = d, h = h)
#' kde_polysph(x = X[1, , drop = FALSE], X = X[-1, ], d = d, h = h, log = TRUE)
#' @export
log_cv_kde_polysph <- function(X, d, h, weights = as.numeric( c()), wrt_unif = FALSE, normalized = TRUE, intrinsic = FALSE, norm_X = FALSE, kernel = 1L, kernel_type = 1L, k = 10.0) {
    .Call(`_polykde_log_cv_kde_polysph`, X, d, h, weights, wrt_unif, normalized, intrinsic, norm_X, kernel, kernel_type, k)
}

#' @title Computation of the softplus function
#'
#' @description Computes the softplus function \eqn{\log(1+e^{t})} in a
#' direct way, as it is supposed to be evaluated for large negative \eqn{t}.
#'
#' @param t vector or matrix.
#' @return The softplus function evaluated at \code{t}.
#' @examples
#' curve(log(polykde:::sfp(rbind(5 * (1 - x)))), from = -10, to = 1)
#' @noRd
sfp <- function(t) {
    .Call(`_polykde_sfp`, t)
}

#' @title Projection onto the polysphere
#'
#' @description Projects points on \eqn{\mathbb{R}^{d_1 + \cdots + d_r + r}}
#' onto the polysphere \eqn{\mathcal{S}^{d_1} \times \cdots \times
#' \mathcal{S}^{d_r}} by normalizing each block of \eqn{d_j} coordinates.
#'
#' @param x a matrix of size \code{c(n, sum(d) + r)}.
#' @param ind_dj \code{0}-based index separating the blocks of spheres that
#' is computed with \code{\link{comp_ind_dj}}.
#' @return A matrix of size \code{c(n, sum(d) + r)} with the projected points.
#' @examples
#' # Example on (S^1)^2
#' d <- c(1, 1)
#' x <- rbind(c(2, 0, 1, 1))
#' proj_polysph(x, ind_dj = comp_ind_dj(d))
#' @export
proj_polysph <- function(x, ind_dj) {
    .Call(`_polykde_proj_polysph`, x, ind_dj)
}

#' @title Polyspherical distance
#'
#' @description Computation of the distance between points \eqn{\boldsymbol{x}}
#' and \eqn{\boldsymbol{y}} on the polysphere
#' \eqn{\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}}:
#' \deqn{\sqrt{\sum_{j=1}^r
#' d_{\mathcal{S}^{d_j}}(\boldsymbol{x}_j, \boldsymbol{y}_j)^2},}
#' where \eqn{d_{\mathcal{S}^{d_j}}(\boldsymbol{x}_j, \boldsymbol{y}_j)=
#' \cos^{-1}(\boldsymbol{x}_j' \boldsymbol{y}_j)}.
#'
#' @inheritParams proj_polysph
#' @param y either a matrix of size \code{c(m, sum(d) + r)} or a vector of
#' length \code{sum(d) + r}.
#' @inheritParams proj_polysph
#' @param norm_x,norm_y ensure a normalization of the data? Default to
#' \code{FALSE}.
#' @param std standardize distance to \eqn{[0,1]}? Uses that the maximum
#' distance is \eqn{\sqrt{r}\pi}. Defaults to \code{TRUE}.
#' @return
#' \itemize{
#' \item{\code{dist_polysph}: a vector of size \code{n} with the distances
#' between \code{x} and \code{y}.}
#' \item{\code{dist_polysph_matrix}: a matrix of size \code{c(n, n)} with the
#' pairwise distances of \code{x}.}
#' \item{\code{dist_polysph_cross}: a matrix of distances of size
#' \code{c(n, m)} with the cross distances between \code{x} and \code{y}.}
#' }
#' @examples
#' # Example on S^2 x S^3 x S^1
#' d <- c(2, 3, 1)
#' ind_dj <- comp_ind_dj(d)
#' n <- 3
#' x <- r_unif_polysph(n = n, d = d)
#' y <- r_unif_polysph(n = n, d = d)
#'
#' # Distances of x to y
#' dist_polysph(x = x, y = y, ind_dj = ind_dj, std = FALSE)
#' dist_polysph(x = x, y = y[1, , drop = FALSE], ind_dj = ind_dj, std = FALSE)
#'
#' # Pairwise distance matrix of x
#' dist_polysph_matrix(x = x, ind_dj = ind_dj, std = FALSE)
#'
#' # Cross distances between x and y
#' dist_polysph_cross(x = x, y = y, ind_dj = ind_dj, std = FALSE)
#' @export
dist_polysph <- function(x, y, ind_dj, norm_x = FALSE, norm_y = FALSE, std = TRUE) {
    .Call(`_polykde_dist_polysph`, x, y, ind_dj, norm_x, norm_y, std)
}

#' @rdname dist_polysph
#' @export
dist_polysph_cross <- function(x, y, ind_dj, norm_x = FALSE, norm_y = FALSE, std = TRUE) {
    .Call(`_polykde_dist_polysph_cross`, x, y, ind_dj, norm_x, norm_y, std)
}

#' @title Diamond cross-products
#'
#' @description Given a matrix \eqn{\boldsymbol{X}} whose \eqn{n} rows are on
#' a polysphere \eqn{\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}
#' \subset \mathbb{R}^p,}:
#' \itemize{
#' \item{\code{diamond_crossprod} computes the \eqn{n\times p\times p} cube
#' whose rows are \eqn{\boldsymbol{X}_i \diamond \boldsymbol{X}_i'},
#' \eqn{i = 1, \ldots, n}, and \eqn{\diamond} is a block-by-block product.}
#' \item{\code{diamond_rcrossprod} computes the \eqn{n\times n\times r} cube
#' formed by \eqn{(\boldsymbol{X}_{ik}' \boldsymbol{X}_{jk})_{ijk}},
#' \eqn{k = 1, \ldots, r}.}
#' }
#'
#' @inheritParams kde_polysph
#' @inheritParams proj_polysph
#' @return
#' \itemize{
#' \item{\code{diamond_crossprod}: an array of size \code{c(nrow(X), ncol(X),
#'  ncol(X))}.}
#' \item{\code{diamond_rcrossprod}: an array of size \code{c(nrow(X), nrow(X),
#' r)}.}
#' }
#' @examples
#' d <- c(1, 2)
#' X <- r_unif_polysph(n = 2, d = d)
#' polykde:::diamond_crossprod(X = X, ind_dj = comp_ind_dj(d))
#' polykde:::diamond_rcrossprod(X = X, ind_dj = comp_ind_dj(d))
#' @noRd
diamond_crossprod <- function(X, ind_dj) {
    .Call(`_polykde_diamond_crossprod`, X, ind_dj)
}

#' @name diamond_crossprod
#' @noRd
diamond_rcrossprod <- function(X, ind_dj) {
    .Call(`_polykde_diamond_rcrossprod`, X, ind_dj)
}

#' @title Symmetrize a matrix
#'
#' @description Symmetrizes a matrix \eqn{\boldsymbol{A}} by returning
#' \eqn{(\boldsymbol{A} + \boldsymbol{A}') / 2}.
#'
#' @param A a matrix.
#' @param add return simply the addition
#' \eqn{\boldsymbol{A} + \boldsymbol{A}'}? Defaults to \code{FALSE}
#' @return A symmetric matrix with the same dimensions as \code{A}.
#' @examples
#' polykde:::s(matrix(rnorm(4), nrow = 2, ncol = 2))
#' @noRd
s <- function(A, add = FALSE) {
    .Call(`_polykde_s`, A, add)
}

#' @title Projection matrices \eqn{\boldsymbol{P}} and \eqn{\boldsymbol{A}}
#'
#' @description Computation of the projection matrices \eqn{\boldsymbol{P}}
#' and \eqn{\boldsymbol{A}}. The \eqn{jj}-block of \eqn{\boldsymbol{P}} is
#' \eqn{\boldsymbol{I}_{d_j} - \boldsymbol{x}_j \boldsymbol{x}_j'}. The
#' \eqn{jj}-block of \eqn{\boldsymbol{A}} is \eqn{(\boldsymbol{x}_j'
#' \boldsymbol{v}_j) \boldsymbol{I}_{d_j}}, \eqn{j=1,\ldots,r}.
#'
#' @param x,v row vectors of size \code{sum(d) + r}.
#' @inheritParams proj_polysph
#' @param orth return the orthogonal complement of \eqn{\boldsymbol{P}},
#' \eqn{\boldsymbol{I} - \boldsymbol{P}}? Defaults to \code{FALSE}.
#' @return A list with the matrices \eqn{\boldsymbol{P}} and
#' \eqn{\boldsymbol{A}}. Both matrices have size
#' \code{c(sum(d) + r, sum(d) + r)}.
#' @examples
#' d <- c(1, 2)
#' x <- r_unif_polysph(n = 1, d = d)
#' v <- r_unif_polysph(n = 1, d = d)
#' polykde:::AP(x = x, v = v, ind_dj = comp_ind_dj(d))
#' @noRd
AP <- function(x, v, ind_dj, orth = FALSE) {
    .Call(`_polykde_AP`, x, v, ind_dj, orth)
}