Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
#include <Rcpp.h>
#include <RcppParallel.h>
// [[Rcpp::depends(RcppParallel)]]
// [[Rcpp::plugins(cpp11)]]
using namespace Rcpp;
// [[Rcpp::export]]
int count_baseline( NumericVector x ){
return std::count( x.begin(), x.end(), 3.0 ) ;
}
// [[Rcpp::export]]
int count_na_rapi( NumericVector x ){
return std::count_if( x.begin(), x.end(), R_IsNA ) ;
}
// [[Rcpp::export]]
int count_na_rcpp( NumericVector x ){
return std::count_if( x.begin(), x.end(), NumericVector::is_na ) ;
}
// [[Rcpp::export]]
int count_na_proposed( NumericVector x ){
// interpret the vector as a vector of 64 bits integers
uint64_t* p = reinterpret_cast<uint64_t*>(x.begin()) ;
// get the 64 bit integer for a quiet NA
uint64_t na = *reinterpret_cast<uint64_t*>(&NA_REAL) ;
// the mask to nuke the 13th bit
uint64_t mask = ~( (uint64_t(1) << 51 ) );
// using a lambda for the actual work
return std::count_if(p, p + x.size(), [na,mask](uint64_t y){
return ( y & mask ) == na ;
}) ;
}
// [[Rcpp::export]]
int par_count_baseline( NumericVector x ){
auto count_chunk = [=](const tbb::blocked_range<double*>& r, int init) -> int {
return init + std::count( r.begin(), r.end(), 3.0 );
} ;
return tbb::parallel_reduce(
tbb::blocked_range<double*>(x.begin(), x.end()),
0,
count_chunk,
[]( int x, int y){ return x+y; }
) ;
}
// [[Rcpp::export]]
int par_count_rapi( NumericVector x ){
auto count_chunk = [=](const tbb::blocked_range<double*>& r, int init) -> int {
return init + std::count_if( r.begin(), r.end(), R_IsNA );
} ;
return tbb::parallel_reduce(
tbb::blocked_range<double*>(x.begin(), x.end()),
0,
count_chunk,
[]( int x, int y){ return x+y; }
) ;
}
// [[Rcpp::export]]
int par_count_na_rcpp( NumericVector x ){
auto count_chunk = [=](const tbb::blocked_range<double*>& r, int init) -> int {
return init + std::count_if( r.begin(), r.end(), NumericVector::is_na );
} ;
return tbb::parallel_reduce(
tbb::blocked_range<double*>(x.begin(), x.end()),
0,
count_chunk,
[]( int x, int y){ return x+y; }
) ;
}
// [[Rcpp::export]]
int par_count_na_proposed( NumericVector x ){
// interpret the vector as a vector of 64 bits integers
uint64_t* p = reinterpret_cast<uint64_t*>(x.begin()) ;
// get the 64 bit integer for a quiet NA
uint64_t na = *reinterpret_cast<uint64_t*>(&NA_REAL) ;
// the mask to nuke the 13th bit
uint64_t mask = ~( (uint64_t(1) << 51 ) );
auto count_chunk = [=](const tbb::blocked_range<uint64_t*>& r, int init) -> int {
return init + std::count_if( r.begin(), r.end(), [=](uint64_t y){
return ( y & mask ) == na ;
}) ;
};
return tbb::parallel_reduce(
tbb::blocked_range<uint64_t*>(p, p + x.size()),
0,
count_chunk,
[]( int x, int y){ return x+y; }
) ;
}
/*** R
library(microbenchmark)
bench <- function(n){
x <- rnorm(n)
x[ sample(n, n/10) ] <- NA
microbenchmark(
R = sum(is.na(x)),
count_baseline = count_baseline(x),
count_na_rapi = count_na_rapi(x),
count_na_rcpp = count_na_rcpp(x),
count_na_proposed = count_na_proposed(x),
par_count_baseline = par_count_baseline(x),
par_count_rapi = par_count_rapi(x),
par_count_na_rcpp = par_count_na_rcpp(x),
par_count_na_proposed = par_count_na_proposed(x)
)
}
# bench(1e5)
# bench(1e6)
*/
@romainfrancois

This comment has been minimized.

Copy link
Owner Author

commented Nov 8, 2017

> bench(1e5)
Unit: microseconds
                  expr     min       lq       mean   median       uq       max neval cld
                     R 154.354 180.7435 1163.07663 204.7635 388.2565 84133.718   100   a
        count_baseline  28.492  38.7965   56.18033  44.6670  47.9475  1048.677   100   a
         count_na_rapi 228.332 299.0865  344.50478 323.1225 360.9430  1181.584   100   a
         count_na_rcpp 315.074 374.7255  422.69063 405.0850 444.1935  1345.352   100   a
     count_na_proposed  59.312  67.2630   84.52644  72.3595  77.4710  1030.615   100   a
    par_count_baseline  23.212  32.7935   81.42861  38.0210  53.2545  1788.471   100   a
        par_count_rapi 110.723 136.3875  166.68644 150.1515 174.6015   871.742   100   a
     par_count_na_rcpp 154.514 179.2490  210.85263 194.4715 217.2000  1070.035   100   a
 par_count_na_proposed  39.007  53.4780   82.55541  61.2895  82.3975  1044.195   100   a


> bench(1e6)
Unit: microseconds
                  expr      min        lq      mean    median        uq      max neval      cld
                     R 1875.361 2112.8835 3188.7986 2846.9530 3183.6135 6783.091   100       g 
        count_baseline  515.749  611.5435  674.6513  656.9820  727.3860 1035.686   100  b      
         count_na_rapi 2389.200 2645.3290 2902.3169 2829.9430 3033.8420 4937.296   100      f  
         count_na_rcpp 3120.384 3525.3630 3895.4203 3829.9835 4057.1865 6455.599   100        h
     count_na_proposed  759.485  853.1500  941.1813  940.3625  997.6140 1717.260   100   c     
    par_count_baseline  293.698  363.0530  402.2780  396.3610  433.1535  581.647   100 a       
        par_count_rapi 1031.302 1153.8500 1414.0163 1356.8435 1629.6385 3686.571   100    d    
     par_count_na_rcpp 1462.124 1614.7535 1797.6803 1783.1680 1918.9350 2548.485   100     e   
 par_count_na_proposed  392.947  470.7030  525.8884  525.2660  576.1705  775.663   100 ab     
 
 
> bench(1e7)
Unit: milliseconds
                  expr       min        lq      mean    median        uq        max neval     cld
                     R 28.160836 29.604265 37.978571 39.757212 41.190611 206.602008   100       g
        count_baseline  5.814813  6.223462  6.822613  6.628946  7.098646  10.003304   100  bc    
         count_na_rapi 23.750143 24.493118 26.412370 25.767943 27.779770  42.620657   100      f 
         count_na_rcpp 31.864211 33.586734 35.574865 35.587514 37.377092  42.523575   100       g
     count_na_proposed  7.956405  8.208651  8.915667  8.416500  9.165068  17.072256   100   cd   
    par_count_baseline  2.803966  3.695622  3.969127  3.797114  3.978292  15.018628   100 a      
        par_count_rapi 10.408895 10.644652 11.532624 11.124726 11.883716  17.623947   100    d   
     par_count_na_rcpp 14.732422 14.907101 15.976666 15.529074 16.575572  21.204417   100     e  
 par_count_na_proposed  4.025824  4.180822  4.653866  4.352212  4.757907   9.028662   100 ab 
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.