r - 快速检查 Rcpp 中的缺失值

标签 r performance rcpp

这个问题链接到NA values in Rcpp conditional .

我基本上有一些循环多个(双)元素的 Rcpp 代码。我需要检查每个元素是否有缺失值(我不能使用矢量化)。让我们计算向量中缺失值的数量，就像最小的可重复示例:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
int nb_na(const NumericVector& x) {
  int n = x.size();
  int c = 0;
  for (int i = 0; i < n; i++) if (R_IsNA(x[i])) c++;
  return c;
}

// [[Rcpp::export]]
int nb_na3(const NumericVector& x) {
  int n = x.size();
  int c = 0;
  for (int i = 0; i < n; i++) if (x[i] == 3) c++;
  return c;
}

// [[Rcpp::export]]
LogicalVector na_real(NumericVector x) {
  return x == NA_REAL;
}

然后，在 R 中，我们得到:
> x <- rep(c(1, 2, NA), 1e4)

> x2 <- replace(x, is.na(x), 3)

> microbenchmark::microbenchmark(
+   nb_na(x),
+   nb_na3(x2)
+ )
Unit: microseconds
       expr     min      lq      mean  median       uq      max neval
   nb_na(x) 135.633 135.982 153.08586 139.753 140.3115 1294.928   100
 nb_na3(x2)  22.490  22.908  30.14005  23.188  23.5025  684.026   100

> all.equal(nb_na(x), nb_na3(x2))
[1] TRUE

> na_real(x[1:3])
[1] NA NA NA

如链接问题所述，您不能只检查 x[i] == NA_REAL因为它总是返回一个缺失值。然而，使用 R_IsNA(x[i])用数值(例如 3 )检查相等性要慢得多。

基本上，我想要一个解决方案，我可以在其中检查单个值是否为缺失值。此解决方案应该与使用数值检查相等性一样快。

最佳答案

检查缺失值或任何 NaN 特定变量总是比检查特定值更昂贵。那只是浮点运算。

但是，您的代码仍有改进的余地。我鼓励您使用 NumericVector::is_na而不是 R_IsNA但这主要是化妆品。

然后分支可能很昂贵，即我会替换 if (R_IsNA(x[i])) c++;来自 c += NumericVector::is_na(x[i]) .这给出了这个版本:

// [[Rcpp::export]]
int nb_na4(const NumericVector& x) {
  int n = x.size();
  int c = 0;
  for (int i = 0; i < n; i++) c += NumericVector::is_na(x[i]) ;
  return c;
}

然后迭代 int并访问 x[i]可以使用 std::count_if algorithm 替换.这是存在的理由。导致这个版本:
// [[Rcpp::export]]
int nb_na5(const NumericVector& x) {
  return std::count_if(x.begin(), x.end(), NumericVector::is_na ) ;
}

现在如果性能仍然不够好，你可能想尝试并行化，为此我通常使用 tbb图书馆来自 RcppParallel包裹。
// [[Rcpp::export]]
int nb_na6(const NumericVector& x) {
  return tbb::parallel_reduce( 
    tbb::blocked_range<const double*>(x.begin(), x.end()),
    0, 
    [](const tbb::blocked_range<const double*>& r, int init) -> int {
      return init + std::count_if( r.begin(), r.end(), NumericVector::is_na );
    }, 
    []( int x, int y){ return x+y; }
  ) ;
}

使用此功能进行基准测试:
library(microbenchmark)

bench <- function(n){
  x <- rep(c(1, 2, NA), n)
  microbenchmark(
    nb_na = nb_na(x), 
    nb_na4 = nb_na4(x), 
    nb_na5 = nb_na5(x), 
    nb_na6 = nb_na6(x)
  )
}
bench(1e5)

在我的机器上，我得到:
> bench(1e4)
Unit: microseconds
expr    min      lq      mean  median       uq     max neval  cld
nb_na  84.358 94.6500 107.41957 110.482 118.9580 137.393   100    d
nb_na4 59.984 69.4925  79.42195  82.442  85.9175 106.567   100  b  
nb_na5 65.047 75.2625  85.17134  87.501  93.0315 116.993   100   c 
nb_na6 39.205 51.0785  59.20582  54.457  68.9625  97.225   100 a   

> bench(1e5)
Unit: microseconds
expr     min       lq     mean   median       uq      max neval  cld
nb_na  730.416 732.2660 829.8440 797.4350 872.3335 1410.467   100    d
nb_na4 520.800 521.6215 598.8783 562.7200 657.1755 1059.991   100  b  
nb_na5 578.527 579.3805 664.8795 626.5530 710.5925 1166.365   100   c 
nb_na6 294.486 345.2050 368.6664 353.6945 372.6205  897.552   100 a

另一种方法是完全绕过浮点运算，并假装向量是 long long 的向量。 ，又名 64 位整数，并将这些值与 NA_REAL 的位模式进行比较:
  > devtools::install_github( "ThinkR-open/seven31" )
  > seven31::reveal(NA, NaN, +Inf, -Inf )
  0 11111111111 ( NaN ) 0000000000000000000000000000000000000000011110100010 : NA
  0 11111111111 ( NaN ) 1000000000000000000000000000000000000000000000000000 : NaN
  0 11111111111 ( NaN ) 0000000000000000000000000000000000000000000000000000 : +Inf
  1 11111111111 ( NaN ) 0000000000000000000000000000000000000000000000000000 : -Inf

使用此 hack 的串行解决方案:
// [[Rcpp::export]]
int nb_na7( const NumericVector& x){
  const long long* p = reinterpret_cast<const long long*>(x.begin()) ;
  long long na = *reinterpret_cast<long long*>(&NA_REAL) ;

  return std::count(p, p + x.size(), na ) ;

}

然后是并行版本:
// [[Rcpp::export]]
int nb_na8( const NumericVector& x){
  const long long* p = reinterpret_cast<const long long*>(x.begin()) ;
  long long na = *reinterpret_cast<long long*>(&NA_REAL) ;

  auto count_chunk = [=](const tbb::blocked_range<const long long*>& r, int init) -> int {
    return init + std::count( r.begin(), r.end(), na);
  } ;

  return tbb::parallel_reduce( 
    tbb::blocked_range<const long long*>(p, p + x.size()),
    0, 
    count_chunk, 
    []( int x, int y){ return x+y; }
  ) ;

}

  > bench(1e5)
  Unit: microseconds
     expr     min       lq     mean   median       uq      max neval    cld
    nb_na 730.346 762.5720 839.9479 857.5865 881.8635 1045.048   100      f
   nb_na4 520.946 521.6850 589.0911 578.2825 653.4950  832.449   100    d  
   nb_na5 578.621 579.3245 640.9772 616.8645 701.8125  890.736   100     e 
   nb_na6 291.115 307.4300 340.1626 344.7955 360.7030  484.261   100   c   
   nb_na7 122.156 123.4990 141.1954 132.6385 149.7895  253.988   100  b    
   nb_na8  69.356  86.9980 109.6427 115.2865 126.2775  182.184   100 a     

  > bench(1e6)
  Unit: microseconds
     expr      min        lq      mean    median        uq      max neval  cld
    nb_na 7342.984 7956.3375 10261.583 9227.7450 10869.605 79757.09   100    d
   nb_na4 5286.970 5721.9150  7659.009 6660.2390  9234.646 31141.47   100   c 
   nb_na5 5840.946 6272.7050  7307.055 6883.2430  8205.117 10420.48   100   c 
   nb_na6 2833.378 2895.7160  3891.745 3049.4160  4054.022 18242.26   100  b  
   nb_na7 1661.421 1791.1085  2708.992 1916.6055  2232.720 60827.63   100 ab  
   nb_na8  650.639  869.6685  1289.373  939.0045  1291.025 10223.29   100 a

这假设只有一个位模式来表示 NA .

这是我的整个文件以供引用:
#include <Rcpp.h>
#include <RcppParallel.h>

// [[Rcpp::depends(RcppParallel)]]
// [[Rcpp::plugins(cpp11)]]
using namespace Rcpp;

// [[Rcpp::export]]
int nb_na(const NumericVector& x) {
  int n = x.size();
  int c = 0;
  for (int i = 0; i < n; i++) if (R_IsNA(x[i])) c++;
  return c;
}

// [[Rcpp::export]]
int nb_na4(const NumericVector& x) {
  int n = x.size();
  int c = 0;
  for (int i = 0; i < n; i++) c += NumericVector::is_na(x[i]) ;
  return c;
}

// [[Rcpp::export]]
int nb_na5(const NumericVector& x) {
  return std::count_if(x.begin(), x.end(), NumericVector::is_na ) ;
}

// [[Rcpp::export]]
int nb_na6(const NumericVector& x) {
  return tbb::parallel_reduce( 
    tbb::blocked_range<const double*>(x.begin(), x.end()),
    0, 
    [](const tbb::blocked_range<const double*>& r, int init) -> int {
      return init + std::count_if( r.begin(), r.end(), NumericVector::is_na );
    }, 
    []( int x, int y){ return x+y; }
  ) ;
}

// [[Rcpp::export]]
int nb_na7( const NumericVector& x){
  const long long* p = reinterpret_cast<const long long*>(x.begin()) ;
  long long na = *reinterpret_cast<long long*>(&NA_REAL) ;

  return std::count(p, p + x.size(), na ) ;

}

// [[Rcpp::export]]
int nb_na8( const NumericVector& x){
  const long long* p = reinterpret_cast<const long long*>(x.begin()) ;
  long long na = *reinterpret_cast<long long*>(&NA_REAL) ;

  auto count_chunk = [=](const tbb::blocked_range<const long long*>& r, int init) -> int {
    return init + std::count( r.begin(), r.end(), na);
  } ;

  return tbb::parallel_reduce( 
    tbb::blocked_range<const long long*>(p, p + x.size()),
    0, 
    count_chunk, 
    []( int x, int y){ return x+y; }
  ) ;

}

/*** R
library(microbenchmark)

bench <- function(n){
  x <- rep(c(1, 2, NA), n)
  microbenchmark(
    nb_na = nb_na(x), 
    nb_na4 = nb_na4(x), 
    nb_na5 = nb_na5(x), 
    nb_na6 = nb_na6(x), 
    nb_na7 = nb_na7(x), 
    nb_na8 = nb_na8(x)
  )
}
bench(1e5)
bench(1e6)
*/

关于r - 快速检查 Rcpp 中的缺失值，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/46892399/

上一篇：cypress - Cypress 是否在每次测试前清除浏览器缓存？

下一篇：performance - Erlang 'catch'表达式与try/catch在效率方面的关系

相关文章：

r - R 中的过滤/子集应用于多列

r - 在 Cassandra 中将宽表转换为长格式

r - 如何使用 3 个具有复杂索引进度的 for 循环来加速此计算？

c++ - 来自头文件/源文件的包含文件

performance - R 与 Matlab : Explanation for speed difference for rnorm, qnorm 和 pnorm 函数

r - 合并具有一个共同元素的集合 R

R:根据条件保持行不重复(多列)

performance - gitlab接口(interface)性能慢

Ubuntu ld 找不到 lRcpp

c++ - 用于子集字符串的 Rcpp 函数

©2024 IT工具网  联系我们