`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/