Rcpp中缺失值的快速检测
Rcpp中缺失值的快速检测
提问于 2017-10-23 14:53:17
这个问题与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)检查等式要慢得多。
基本上，我想要一个解决方案，在这里我可以检查单个值是否是一个缺失的值。此解决方案应该与使用数字值检查等式一样快。。
回答 2
Stack Overflow用户
回答已采纳
发布于 2017-10-24 12:01:36
检查缺少的值或任何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]可以用 algorithm代替。这是存在的理由。导致这一版本：
// [[Rcpp::export]]
int nb_na5(const NumericVector& x) {
return std::count_if(x.begin(), x.end(), NumericVector::is_na ) ;
}现在，如果性能仍然不够好，您可能需要尝试并行化，为此，我通常使用来自RcppParallel包的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使用此黑客的串行解决方案：
// [[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)
*/Stack Overflow用户
发布于 2017-10-23 15:08:19
检查(IEEE)缺少的浮点值是一种昂贵的操作，没有办法避免它。这与R.
这就是为什么我们对即将到来的R中的ALTREP感到兴奋的一个原因--例如，我们可以跟踪一个双/实向量是否包含缺失的值--如果没有，我们就不用浪费时间去寻找它们。虽然没有更新以提及ALTREP，但您可以从https://github.com/HenrikBengtsson/Wishlist-for-R/issues/12获得gist
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/46892399
复制相关文章
相似问题
腾讯云开发者
