Transferring an order to Rcpp, i.e. Base :: order ()

Question

I have a ton of code using the base :: order () command, and I am very lazy to encode the code in rcpp. Since Rcpp only supports sorting, but not ordering, I spent 2 minutes creating this function:

// [[Rcpp::export]] Rcpp::NumericVector order_cpp(Rcpp::NumericVector invec){ int leng = invec.size(); NumericVector y = clone(invec); for(int i=0; i<leng; ++i){ y[sum(invec<invec[i])] = i+1; } return(y); } 

It works somehow. If the vectors contain unique numbers, I get the same result as order (). If they are not unique, the results are different, but not erroneous (in fact, there is no unique solution).

Using it:

 c=sample(1:1000,500) all.equal(order(c),order_cpp(c)) microbenchmark(order(c),order_cpp(c)) Unit: microseconds expr min lq median uq max neval order(c) 33.507 36.223 38.035 41.356 78.785 100 order_cpp(c) 2372.889 2427.071 2466.312 2501.932 2746.586 100 

Oh! I need an efficient algorithm. Ok, so I dug up the bubblesort implementation and adapted it:

  // [[Rcpp::export]] Rcpp::NumericVector bubble_order_cpp2(Rcpp::NumericVector vec){ double tmp = 0; int n = vec.size(); Rcpp::NumericVector outvec = clone(vec); for (int i = 0; i <n; ++i){ outvec[i]=static_cast<double>(i)+1.0; } int no_swaps; int passes; passes = 0; while(true) { no_swaps = 0; for (int i = 0; i < n - 1 - passes; ++i) { if(vec[i] > vec[i+1]) { no_swaps++; tmp = vec[i]; vec[i] = vec[i+1]; vec[i+1] = tmp; tmp = outvec[i]; outvec[i] = outvec[i+1]; outvec[i+1] = tmp; }; }; if(no_swaps == 0) break; passes++; }; return(outvec); } 

Well, this is better - but not great:

 microbenchmark(order(c),order_cpp(c),bubble_order_cpp2(c),sort(c),c[order(c)]) Unit: microseconds expr min lq median uq max neval order(c) 33.809 38.034 40.1475 43.3170 72.144 100 order_cpp(c) 2339.080 2435.675 2478.5385 2526.8350 3535.637 100 bubble_order_cpp2(c) 219.752 231.977 234.5430 241.1840 322.383 100 sort(c) 59.467 64.749 68.2205 75.4645 148.815 100 c[order(c)] 38.336 41.204 44.3735 48.1460 93.878 100 

Another conclusion: it is ordered faster than sorted.

Ok, then for shorter vectors:

 c=sample(1:100) microbenchmark(order(c),order_cpp(c),bubble_order_cpp2(c),sort(c),c[order(c)]) Unit: microseconds expr min lq median uq max neval order(c) 10.566 11.4710 12.8300 14.1880 63.089 100 order_cpp(c) 95.689 100.8200 102.7825 107.3105 198.018 100 bubble_order_cpp2(c) 9.962 11.1700 12.0750 13.2830 64.598 100 sort(c) 39.242 41.5065 42.5620 46.3355 155.758 100 c[order(c)] 11.773 12.6790 13.5840 15.9990 82.710 100 

Ok, I missed the RcppArmadillo function:

 // [[Rcpp::export]] Rcpp::NumericVector ordera(arma::vec x) { return(Rcpp::as<Rcpp::NumericVector>(wrap(arma::sort_index( x )+1)) ); } microbenchmark(order(c),order_(c),ordera(c)) Unit: microseconds expr min lq median uq max neval order(c) 9.660 11.169 11.773 12.377 46.185 100 order_(c) 4.529 5.133 5.736 6.038 34.413 100 ordera(c) 4.227 4.830 5.434 6.038 60.976 100 

Answer 1

Here is a simple version that uses Rcpp sugar to implement the order function. We check for duplicates to ensure that everything works as expected. (There is also an error with the Rcpp sort method when there are NA s, so you might also want to check if this will be fixed by the next version).

 #include <Rcpp.h> using namespace Rcpp; // [[Rcpp::export]] IntegerVector order_(NumericVector x) { if (is_true(any(duplicated(x)))) { Rf_warning("There are duplicates in 'x'; order not guaranteed to match that of R base::order"); } NumericVector sorted = clone(x).sort(); return match(sorted, x); } /*** R library(microbenchmark) x <- runif(1E5) identical( order(x), order_(x) ) microbenchmark( order(x), order_(x) ) */ 

gives me

 > Rcpp::sourceCpp('~/test-order.cpp') > set.seed(456) > library(microbenchmark) > x <- runif(1E5) > identical( order(x), order_(x) ) [1] TRUE > microbenchmark( + order(x), + order_(x) + ) Unit: milliseconds expr min lq median uq max neval order(x) 15.48007 15.69709 15.86823 16.21142 17.22293 100 order_(x) 10.81169 11.07167 11.40678 11.87135 48.66372 100 > 

Of course, if you are comfortable with an output that does not match R, you can remove the duplicate check - x[order_(x)] will still be sorted correctly; more specifically, all(x[order(x)] == x[order_(x)]) should return TRUE .

Answer 2

Another solution based on C ++ 11:

 // [[Rcpp::plugins(cpp11)]] #include <Rcpp.h> using namespace Rcpp; template <int RTYPE> IntegerVector order_impl(const Vector<RTYPE>& x, bool desc) { auto n = x.size(); IntegerVector idx = no_init(n); std::iota(idx.begin(), idx.end(), static_cast<size_t>(1)); if (desc) { auto comparator = [&x](size_t a, size_t b){ return x[a - 1] > x[b - 1]; }; std::stable_sort(idx.begin(), idx.end(), comparator); } else { auto comparator = [&x](size_t a, size_t b){ return x[a - 1] < x[b - 1]; }; std::stable_sort(idx.begin(), idx.end(), comparator); // simulate na.last size_t nas = 0; for (size_t i = 0; i < n; ++i, ++nas) if (!Vector<RTYPE>::is_na(x[idx[i] - 1])) break; std::rotate(idx.begin(), idx.begin() + nas, idx.end()); } return idx; } // [[Rcpp::export]] IntegerVector order2(SEXP x, bool desc = false) { switch(TYPEOF(x)) { case INTSXP: return order_impl<INTSXP>(x, desc); case REALSXP: return order_impl<REALSXP>(x, desc); case STRSXP: return order_impl<STRSXP>(x, desc); default: stop("Unsupported type."); } } /***R int <- sample.int(1000, 1E5, replace = TRUE) dbl <- runif(1E5) chr <- sample(letters, 1E5, replace = TRUE) library(benchr) benchmark(order(int), order2(int)) benchmark(order(dbl), order2(dbl)) benchmark(order(chr), order2(chr)) */ 

Compare performance:

 R> int <- sample.int(1000, 1E5, replace = TRUE) R> dbl <- runif(1E5) R> chr <- sample(letters, 1E5, replace = TRUE) R> library(benchr) R> benchmark(order(int), order2(int)) Benchmark summary: Time units : microseconds expr n.eval min lw.qu median mean up.qu max total relative order(int) 100 442 452 464 482 486 1530 48200 1.0 order2(int) 100 5150 5170 5220 5260 5270 6490 526000 11.2 R> benchmark(order(dbl), order2(dbl)) Benchmark summary: Time units : milliseconds expr n.eval min lw.qu median mean up.qu max total relative order(dbl) 100 13.90 14.00 14.20 14.80 15.8 17.4 1480 1.98 order2(dbl) 100 7.11 7.13 7.15 7.26 7.3 8.8 726 1.00 R> benchmark(order(chr), order2(chr)) Benchmark summary: Time units : milliseconds expr n.eval min lw.qu median mean up.qu max total relative order(chr) 100 128.0 131.0 133.0 133.0 134.0 148.0 13300 7.34 order2(chr) 100 17.7 17.9 18.1 18.2 18.3 22.2 1820 1.00 

Note that the radix method from the order base is much faster.

Answer 3

Here is another approach using std::sort .

 typedef std::pair<int, double> paired; bool cmp_second(const paired & left, const paired & right) { return left.second < right.second; } Rcpp::IntegerVector order(const Rcpp::NumericVector & x) { const size_t n = x.size(); std::vector<paired> pairs; pairs.reserve(n); for(size_t i = 0; i < n; i++) pairs.push_back(std::make_pair(i, x(i))); std::sort(pairs.begin(), pairs.end(), cmp_second<paired>); Rcpp::IntegerVector result = Rcpp::no_init(n); for(size_t i = 0; i < n; i++) result(i) = pairs[i].first; return result; }