nacount.cpp

#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)
*/

> 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