`r knitr::opts_chunk\$set(tidy=FALSE)` # Extended example: Rcpp Indices of the largest `n` elements of a large vector (from [StackOverflow][]). # Brainstorm: possible solutions? **Ideas??** # Implementations ## R solutions Sort in decreasing order, finding the nth largest, then scan for values >= the largest value. ```{r r-solutions-1} f0 <- function(x, n) which(x >= sort(x, decreasing=TRUE)[n], arr.ind=TRUE) ``` Partially sort, accomodating the limitation that partial sorting is only supported in increasing order ```{r r-solutions-2} f1 <- function(x, n) which(x >= -sort(-x, partial=n)[n], arr.ind=TRUE) ``` ## An Rcpp / Standard Template Library solution Start by including the appropriate headers and using namespaces for convenience #include #include using namespace Rcpp; using namespace std; Arrange to expose the C++ function to R // [[Rcpp::export]] IntegerVector top_i_pq(NumericVector v, int n) and define some variables, most importantly a `priority_queue` to hold as a pair the numeric value and index. The queue is ordered so the smallest values are at the 'top', with small relying on the standard `pair<>` comparator. typedef pair Elt; priority_queue< Elt, vector, greater > pq; vector result; Walk through the input data, adding it to the queue if either (a) there are fewere than `n` values in the queue or (b) the current value is larger than the smallest value in the queue. In the latter case, pop off the smallest value, and insert it's replacement. In this way the priority queue always contains the `n` largest elements. for (int i = 0; i != v.size(); ++i) { if (pq.size() < n) pq.push(Elt(v[i], i)); else { Elt elt = Elt(v[i], i); if (pq.top() < elt) { pq.pop(); pq.push(elt); } } } Finally pop the indexes from the priority queue into the return vector, remembering to translate to 1-based R coordinates. result.reserve(pq.size()); while (!pq.empty()) { result.push_back(pq.top().second + 1); pq.pop(); } and return the result to R return wrap(result); The full code is available at ```{r top_i_pq-source, eval=FALSE} system.file(package="UnderstandingRBioc", "extdata", "top_i_pq.cpp") ``` This has nice memory use (the priority queue and return vector are both small relative to the original data) and should be fast. Problems with this code include: 1. The default comparator `greater` works so that, in the case of a tie spanning the value of the _n_th element, the last, rather than first, duplicate is retained. 2. NA values (and non-finite values?) may not be handled correctly; I'm not sure whether this is true or not (they are handled correctly, though perhaps more by accident than by design). 3. The function only works for `NumericVector` input, but the logic is appropriate for any R data type for which an appropriate ordering relationship is defined. These issues have been taken up, with a faster and more general algorithm, in an [Rcpp Gallery post][RcppGallery]. ## Output and performance ```{r final-output} library(microbenchmark) library(Rcpp) fl <- system.file(package="UnderstandingRBioc", "extdata", "top_i_pq.cpp") sourceCpp(fl) n <- 10000 x <- rnorm(1e7) res0 <- f0(x, n); res1 <- f1(x, n); res2 <- top_i_pq(x, n) identical(res0, res1) identical(res0, res2) identical(res0, sort(res2)) times <- microbenchmark(f0(x, n), f1(x, n), top_i_pq(x, n), times=5) print(times) plot(times) ``` ## Packaging Shared libraries are often most useful when incorporated into R packages. Packages make it easy to distribute shared libraries for re-use by colleagues, to document functionality in a way that is accessible to non-C++ programmers, and to test the code across platforms. Create the skeleton of an Rcpp package with Rcpp.package.skeleton("CppTest") This produces the basic structure of an Rcpp package in a directory `CppTest`; the details need to be completed as outlined in the file `CppTest/Read-and-delete-me`. The `DESCRIPTION` file includes the line Depends: Rcpp (>= 0.10.4) so that the Rcpp package (especially headers) are available when the package is built and installed. The `NAMESPACE` file indicates that a shared library is present and will be used: useDynLib(TestCpp) The `src/Makevars` and `src/Makevars.win` (for Windows) files. The `Makevars` files configure the compiler to find Rcpp headers: PKG_LIBS = `\$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"` The `src/rcpp_hello_world.h` contains the line RcppExport SEXP rcpp_hello_world(); to indicate that the function `rcpp_hello_world` will be visible from within R. If linking to additional libraries is required, add lines like PKG_CPPFLAGS=-I/opt/boost_1_54_0/include -I../inst/include PKG_LIBS = -L/opt/boost_1_54_0/lib -lboost_regex \ `\$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"` to `Makevars`. [StackOverflow]: http://stackoverflow.com/questions/14821859 [RcppGallery]: http://gallery.rcpp.org/articles/top-elements-from-vectors-using-priority-queue/