R/RcppExports.R

In quitefastmst: Euclidean and Mutual Reachability Minimum Spanning Trees

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

#' @title
#' Get or Set the Number of Threads
#'
#' @description
#' These functions get or set the maximal number of OpenMP threads that
#' can be used by \code{\link{knn_euclid}} and \code{\link{mst_euclid}},
#' amongst others.
#'
#' @param n_threads maximal number of threads to use
#'
#' @return
#' \code{omp_get_max_threads} returns the maximal number
#' of threads that will be used during the next call to a parallelised
#' function, not the maximal number of threads possibly available.
#' It there is no built-in support for OpenMP, 1 is always returned.
#'
#' For \code{omp_set_num_threads}, the previous value of \code{max_threads}
#' is returned.
#'
#'
#' @rdname omp
#' @encoding UTF-8
#' @export
omp_set_num_threads <- function(n_threads) {
    .Call(`_quitefastmst_Romp_set_num_threads`, n_threads)
}

#' @rdname omp
#' @export
omp_get_max_threads <- function() {
    .Call(`_quitefastmst_Romp_get_max_threads`)
}

#' @title Euclidean Nearest Neighbours
#'
#' @description
#' If \code{Y} is \code{NULL}, then the function determines the first \code{k}
#' nearest neighbours of each point in \code{X} with respect
#' to the Euclidean distance. It is assumed that each query point is
#' not its own neighbour.
#'
#' Otherwise, for each point in \code{Y}, this function determines the \code{k}
#' nearest points thereto from \code{X}.
#'
#' @details
#' The implemented algorithms, see the \code{algorithm} parameter, assume
#' that \eqn{k} is rather small.
#'
#' Our implementation of K-d trees (Bentley, 1975) has been quite optimised;
#' amongst others, it has good locality of reference (at the cost of making
#' a copy of the input dataset), features the sliding
#' midpoint (midrange) rule suggested by Maneewongvatana and Mound (1999),
#' node pruning strategies inspired by some ideas from (Sample et al., 2001),
#' and a couple of further tuneups proposed by the current author.
#' Still, it is well-known that K-d trees perform well only in spaces of low
#' intrinsic dimensionality.  Thus, due to the so-called curse of
#' dimensionality, for high \code{d}, the brute-force algorithm is recommended.
#'
#' The number of threads is controlled via the \code{OMP_NUM_THREADS}
#' environment variable or via the \code{\link{omp_set_num_threads}} function
#' at runtime. For best speed, consider building the package
#' from sources using, e.g., \code{-O3 -march=native} compiler flags.
#'
#'
#' @references
#' J.L. Bentley, Multidimensional binary search trees used for associative
#' searching, \emph{Communications of the ACM} 18(9), 509–517, 1975,
#' \doi{10.1145/361002.361007}
#'
#' S. Maneewongvatana, D.M. Mount, It's okay to be skinny, if your friends
#' are fat, \emph{4th CGC Workshop on Computational Geometry}, 1999
#'
#' N. Sample, M. Haines, M. Arnold, T. Purcell, Optimizing search
#' strategies in K-d Trees, \emph{5th WSES/IEEE Conf. on Circuits, Systems,
#' Communications & Computers} (CSCC'01), 2001
#'
#'
#' @param X the "database"; a matrix of shape \eqn{n\times d}
#' @param k requested number of nearest neighbours
#' @param Y the "query points"; \code{NULL} or a matrix of shape \eqn{m\times d};
#'     note that setting \code{Y=X}, contrary to \code{NULL},
#'     will include the query points themselves amongst their own neighbours
#' @param algorithm \code{"auto"}, \code{"kd_tree"} or \code{"brute"};
#'     K-d trees can be used for \code{d} between 2 and 20 only;
#'     \code{"auto"} selects \code{"kd_tree"} in low-dimensional spaces
#' @param max_leaf_size maximal number of points in the K-d tree leaves;
#'     smaller leaves use more memory, yet are not necessarily faster;
#'     use \code{0} to select the default value, currently set to 32
#' @param squared whether the output \code{nn.dist} should be based on
#'     the squared Euclidean distance
#' @param verbose whether to print diagnostic messages
#'
#'
#' @return
#' A list with two elements, \code{nn.index} and \code{nn.dist}, is returned.
#'
#' \code{nn.dist} and \code{nn.index} have shape \eqn{n\times k}
#' or \eqn{m\times k}, depending whether \code{Y} is given.
#'
#' \code{nn.index[i,j]} is the index (between \eqn{1} and \eqn{n})
#' of the \eqn{j}-th nearest neighbour of \eqn{i}.
#'
#' \code{nn.dist[i,j]} gives the weight of the edge \code{{i, nn.index[i,j]}},
#' i.e., the distance between the \eqn{i}-th point and its \eqn{j}-th
#' nearest neighbour, \eqn{j=1,\dots,k}.
#' \code{nn.dist[i,]} is sorted nondecreasingly for all \eqn{i}.
#'
#'
#'
#' @examples
#' library("datasets")
#' data("iris")
#' X <- jitter(as.matrix(iris[1:2]))  # some data
#' neighbours <- knn_euclid(X, 1)  # 1-NNs of each point
#' plot(X, asp=1, las=1)
#' segments(X[,1], X[,2], X[neighbours$nn.index,1], X[neighbours$nn.index,2])
#'
#' knn_euclid(X, 5, matrix(c(6, 4), nrow=1))  # five closest points to (6, 4)
#'
#'
#' @seealso \code{\link{mst_euclid}}
#'
#' @rdname knn_euclid
#' @encoding UTF-8
#' @export
knn_euclid <- function(X, k = 1L, Y = NULL, algorithm = "auto", max_leaf_size = 0L, squared = FALSE, verbose = FALSE) {
    .Call(`_quitefastmst_knn_euclid`, X, k, Y, algorithm, max_leaf_size, squared, verbose)
}

#' @title Euclidean and Mutual Reachability Minimum Spanning Trees
#'
#' @description
#' The function determines the/a(*) minimum spanning tree (MST) of a set
#' of \eqn{n} points, i.e., an acyclic undirected connected graph whose
#' vertices represent the points, and edges are weighted by the distances
#' between point pairs and have minimal total weight.
#'
#' MSTs have many uses in, amongst others, topological data analysis
#' (clustering, density estimation, dimensionality reduction,
#' outlier detection, etc.).
#'
#' In clustering and density estimation, the parameter \code{M} plays the role
#' of a smoothing factor; for discussion, see (Campello et al., 2015)
#' and the references therein.
#'
#' For \eqn{M\leq 1}, we get a spanning tree that minimises the sum of
#' Euclidean distances between the points, i.e., the classic Euclidean minimum
#' spanning tree (EMST).   If \eqn{M=1}, the function additionally returns
#' the distance to each point's nearest neighbour.
#'
#' If \eqn{M>1}, the spanning tree is the smallest with respect to
#' the degree-\eqn{M} mutual reachability distance (Campello et al., 2013) given by
#' \eqn{d_M(i, j)=\max\{ c_M(i), c_M(j), d(i, j)\}}, where \eqn{d(i,j)}
#' is the standard Euclidean distance between the \eqn{i}-th and the \eqn{j}-th point,
#' and \eqn{c_M(i)} is the \eqn{i}-th \eqn{M}-core distance defined as the distance
#' between the \eqn{i}-th point and its \eqn{M}-th nearest neighbour
#' (not including the query point itself).
#'
#' Note that (Campello et al., 2013) defines the core distance as the
#' distance to the \eqn{(M-1)}-th nearest neighbour (or the \eqn{M}-th one,
#' but including self).
#'
#'
#' @details
#' (*) Note that if there are many pairs of equidistant points,
#' there can be many minimum spanning trees. In particular, it is likely
#' that there are point pairs with the same mutual reachability distances.
#'
#' To make the definition unambiguous, the \code{mutreach_ties} argument
#' indicates the preference towards connecting to farther/closer points with
#' respect to the original metric, or having smaller/larger core distances
#' in cases of tied distances; see (Gagolewski, 2026).  Empirically,
#' \code{mutreach_ties="dcore_min"} and \code{mutreach_leaves="reconnect_dcore_min"}
#' leads to MSTs with more leaves and hubs.  This is only available
#' in the brute force method.
#'
#' The implemented algorithms, see the \code{algorithm} parameter, assume
#' that \eqn{M} is rather small.
#'
#' Our implementation of K-d trees (Bentley, 1975) has been quite optimised;
#' amongst others, it has good locality of reference (at the cost of making
#' a copy of the input dataset), features the sliding
#' midpoint (midrange) rule suggested by Maneewongvatana and Mound (1999),
#' node pruning strategies inspired by some ideas from (Sample et al., 2001),
#' and a couple of further tuneups proposed by the current author.
#'
#' The "single-tree" version of the Borůvka algorithm is parallelised:
#' in every iteration, it seeks each point's nearest "alien",
#' i.e., the nearest point thereto from another cluster.
#' The "dual-tree" Borůvka version of the algorithm is, in principle, based
#' on (March et al., 2010). As far as our implementation is concerned,
#' the dual-tree approach is often only faster in 2- and 3-dimensional spaces,
#' for \eqn{M\leq 1}, and in a single-threaded setting.  For another
#' (approximate) adaptation of the dual-tree algorithm to mutual
#' reachability distances, see (McInnes and Healy, 2017).
#'
#' The "sesqui-tree" variant (by the current author) is a mixture of the two
#' approaches:  it compares leaves against the full tree and can be run
#' in parallel.  It is usually faster than the single- and dual-tree methods
#' in very low dimensional spaces and usually not much slower than
#' the single-tree variant otherwise.
#'
#' Nevertheless, it is well-known that K-d trees perform well only in spaces
#' of low intrinsic dimensionality (the "curse").  For high \eqn{d},
#' the "brute-force" algorithm is recommended.  Here, we provided a
#' parallelised (see Olson, 1995) version of the Jarník (1930) (a.k.a.
#' Prim, 1957) algorithm, where the distances are computed
#' on the fly (only once for \eqn{M\leq 1}).
#'
#' The number of threads used is controlled via the \code{OMP_NUM_THREADS}
#' environment variable or via the \code{\link{omp_set_num_threads}} function
#' at runtime. For best speed, consider building the package
#' from sources using, e.g., \code{-O3 -march=native} compiler flags.
#'
#'
#' @references
#' V. Jarník, O jistém problému minimálním,
#' \emph{Práce Moravské Přírodovědecké Společnosti} 6, 1930, 57–63.
#'
#' C.F. Olson, Parallel algorithms for hierarchical clustering,
#' Parallel Computing 21(8), 1995, 1313–1325.
#'
#' R. Prim, Shortest connection networks and some generalizations,
#' \emph{The Bell System Technical Journal} 36(6), 1957, 1389–1401.
#'
#' O. Borůvka, O jistém problému minimálním, \emph{Práce Moravské
#' Přírodovědecké Společnosti} 3, 1926, 37–58.
#'
#' W.B. March, R. Parikshit, A.G. Gray, Fast Euclidean minimum spanning
#' tree: Algorithm, analysis, and applications, \emph{Proc. 16th ACM SIGKDD
#' Intl. Conf. Knowledge Discovery and Data Mining (KDD '10)}, 2010, 603–612.
#'
#' J.L. Bentley, Multidimensional binary search trees used for associative
#' searching, \emph{Communications of the ACM} 18(9), 509–517, 1975,
#' \doi{10.1145/361002.361007}
#'
#' S. Maneewongvatana, D.M. Mount, It's okay to be skinny, if your friends
#' are fat, \emph{4th CGC Workshop on Computational Geometry}, 1999
#'
#' N. Sample, M. Haines, M. Arnold, T. Purcell, Optimizing search
#' strategies in K-d Trees, \emph{5th WSES/IEEE Conf. on Circuits, Systems,
#' Communications & Computers} (CSCC'01), 2001
#'
#' R.J.G.B. Campello, D. Moulavi, J. Sander, Density-based clustering based
#' on hierarchical density estimates, \emph{Lecture Notes in Computer Science}
#' 7819, 2013, 160–172. \doi{10.1007/978-3-642-37456-2_14}
#'
#' R.J.G.B. Campello, D. Moulavi, A. Zimek, J. Sander, Hierarchical
#' density estimates for data clustering, visualization, and outlier detection,
#' \emph{ACM Transactions on Knowledge Discovery from Data (TKDD)} 10(1),
#' 2015, 1–51, \doi{10.1145/2733381}
#'
#' L. McInnes, J. Healy, Accelerated hierarchical density-based
#' clustering, \emph{IEEE Intl. Conf. Data Mining Workshops (ICMDW)}, 2017,
#' 33–42, \doi{10.1109/ICDMW.2017.12}
#'
#' M. Gagolewski, quitefastmst, in preparation, 2026, TODO
#'
#'
#' @param X the "database"; a matrix of shape \eqn{n\times d}
#' @param M the smoothing factor a.k.a. the degree of the mutual reachability
#'    distance; \eqn{M\leq 1} gives the ordinary Euclidean distance
#' @param algorithm \code{"auto"}, \code{"single_kd_tree"},
#'    \code{"sesqui_kd_tree"}, \code{"dual_kd_tree"}, or \code{"brute"};
#'    K-d trees can only be used for \eqn{d} between 2 and 20 only;
#'    \code{"auto"} selects \code{"sesqui_kd_tree"} for \eqn{d\leq 20}.
#'    \code{"brute"} is used otherwise
#' @param max_leaf_size maximal number of points in the K-d tree leaves;
#'    smaller leaves use more memory, yet are not necessarily faster;
#'    use \code{0} to select the default value, currently set to 32 for the
#'    single-tree and sesqui-tree and 8 for the dual-tree Borůvka algorithm
#' @param first_pass_max_brute_size minimal number of points in a node to
#'    treat it as a leaf (unless it actually is a leaf) in the first
#'    iteration of the algorithm; use \code{0} to select the default value,
#'    currently set to 32
#' @param mutreach_ties adjustment for mutual reachability distance ambiguity
#'    (for \eqn{M>1} and algorithm \code{"brute"}); one of \code{"dcore_min"},
#'    \code{"dist_max"}, \code{"dist_min"} (default), or \code{"dcore_max"}
#' @param mutreach_leaves a way to postprocess the leaves of the computed tree;
#'    one of \code{"keep"} (default: do nothing),
#'    or \code{"reconnect_dcore_min"} (try reconnecting leaves to
#'    inner vertices which have them amongst their M nearest neighbours;
#'    prefer vertices of the smallest core distance)
#' @param verbose whether to print diagnostic messages
#'
#'
#' @return
#' A list with two $(M=0)$ or four $(M>0)$ elements, \code{mst.index} and
#' \code{mst.dist}, and additionally \code{nn.index} and \code{nn.dist}.
#'
#' \code{mst.index} is a matrix with \eqn{n-1} rows and \eqn{2} columns,
#' whose rows define the tree edges.
#'
#' \code{mst.dist} is a vector of length
#' \eqn{n-1} giving the weights of the corresponding edges.
#'
#' The tree edges are ordered with respect to weights nondecreasingly, and then by
#' the indexes (lexicographic ordering of the \code{(weight, index1, index2)}
#' triples).  For each \code{i}, it holds \code{mst_ind[i,1]<mst_ind[i,2]}.
#'
#' \code{nn.index} is an \eqn{n} by \eqn{M} matrix giving the indexes
#' of each point's nearest neighbours with respect to the Euclidean distance.
#' \code{nn.dist} provides the corresponding Euclidean distances.
#'
#'
#' @examples
#' library("datasets")
#' data("iris")
#' X <- jitter(as.matrix(iris[1:2]))  # some data
#' T <- mst_euclid(X)                 # Euclidean MST of X
#' plot(X, asp=1, las=1)
#' segments(X[T$mst.index[, 1], 1], X[T$mst.index[, 1], 2],
#'          X[T$mst.index[, 2], 1], X[T$mst.index[, 2], 2])
#'
#' @seealso \code{\link{knn_euclid}}
#'
#' @rdname mst_euclid
#' @encoding UTF-8
#' @export
mst_euclid <- function(X, M = 0L, algorithm = "auto", max_leaf_size = 0L, first_pass_max_brute_size = 0L, mutreach_ties = "dist_min", mutreach_leaves = "keep", verbose = FALSE) {
    .Call(`_quitefastmst_mst_euclid`, X, M, algorithm, max_leaf_size, first_pass_max_brute_size, mutreach_ties, mutreach_leaves, verbose)
}