R：大向量的高效迭代子集和过滤

我想更快地执行以下操作.

我有一个由4个元素1, 2, 3, 4组成的向量.我也有一个相同长度的阈值向量1.1, 3.1, 4.1, 5.1.我希望每个元素都找到前next个元素的索引高于相应的阈值.在本例中，我的预期输出为

2, 3, NA, NA:

• 高于阈值1.1的第一个元素after the first one (included)处于索引2(值2).
• 第二个阈值3.1以上的第一个元素的值为4，是第third个元素after the current one at index 2(包括在内).

Base implementation个

start <- Sys.time()
bigg <- rnorm(25000)
thresh <- bigg+0.5
result <- rep(NA, length(bigg))
for(i in 1:length(bigg)) {
  result[i] <- which(bigg[(i+1):length(bigg)]>thresh[i])[1] # the first next element that is higher than thresh
  if(i%%1000==0) print(paste0(i, " ", round(i/length(bigg),3)))
}
end <- Sys.time()
end-start
head(result)

基本上，取满足阈值条件的向量xafter the current one的第一个元素.

我试着用Rcpp

// [[Rcpp::export]]
int cppnextup_(NumericVector x, double thresh, bool is_up = true) {
  int n = x.size();
  //int idx = 0;
  int res = -1;
  for(int idx = 0; idx < n; ++idx) {
    if(x[idx]>thresh && is_up == true) {
      res = idx;
//Rcout << "The value of idx : " << idx <<" "<< x[idx]<<"\n";
      break;
    }
    if(x[idx]<thresh && is_up == false) {
      res = idx;
      //Rcout << "The value of idx : " << idx <<" "<< x[idx]<<"\n";
      break;
    }
  }
  return res;
}

Benchmarking:个

# base --------------------------------------------------------------------

base_ <- function() {
  for(i in 1:length(bigg)) {
    result[i] <- which(bigg[(i+1):length(bigg)]>thresh[i])[1] # the first next element that is higher than thresh
    if(i%%1000==0) print(paste0(i, " ", round(i/length(bigg),3)))
  }
}

# cpp ----------------------------------------------------------------

result_cpp <- rep(NA, length(bigg))
cpp_ <- function() {
  for(i in 1:length(bigg)) {
    result_cpp[i] <- cppnextup_(bigg[(i+1):length(bigg)], thresh[i]) # the first next element that is higher than thresh
    if(i%%1000==0) print(paste0(i, " ", round(i/length(bigg),3)))
  }
}

#result_cpp <- ifelse(result_cpp==-1, NA, result_cpp)
#result_cpp <- result_cpp+1
#all.equal(result, result_cpp)
#[1] TRUE

# benchmark ---------------------------------------------------------------

microbenchmark::microbenchmark(base_(),
                               cpp_(), times=3)
Unit: milliseconds
    expr      min        lq      mean    median        uq       max neval
 base_() 2023.510 2030.3154 2078.7867 2037.1211 2106.4252 2175.7293     3
  cpp_()  661.277  665.3456  718.8851  669.4141  747.6891  825.9641     3

我的Rcpp实现将基本时间减少了65%，有没有更好的(矢量化)方法？寻找任何后端，可以是Rcpp,data.table,dtplyr等.

我的dtplyr次try 得到了所有的NA分:

library(dtplyr)
nx <- length(bigg)
df <- tibble(bigg, thresh)
bigg %>% lazy_dt() %>% mutate(res = which(bigg[row_number():nx]>thresh)[1])
Warning message:
In seq_len(.N):..nx :
  numerical expression has 25000 elements: only the first used

干杯

顺便说一下，我的真实向量有8,406,600个元素.

编辑:矢量化Rcpp

我还有另一个更快的Rcpp函数，它依赖于第一个函数:

// [[Rcpp::export]]
NumericVector cppnextup(NumericVector x, double threshup, bool is_up = true) {
  int n = x.size();
  NumericVector up(n);
  if(is_up == true) {
    up = x + threshup;
  } else {
    up = x - threshup;
  }
//  Rcout << "The value of up : " << up[0] <<" "<< up[1] <<"\n";
  NumericVector result(n);
  int idx = 0;
  for(int i = 0; i < n; ++i) {
    double thisup = up[idx];
    NumericVector thisvect = x[Rcpp::Range((idx), (n-1))];
    
//Rcout <<idx<< " " << "thisvect : " << thisvect[0] <<" thisup: "<< thisup <<" buy " << buy << "\n";
    
    int resi = cppnextup_(thisvect, thisup, is_up = is_up);
    if(resi != 0) {
      result[idx] = resi+1;
    } else {
      result[idx] = resi;
    }
    
    //Rcout << "RESI: " << resi <<" "<< up[1] <<"\n";
    idx = idx + 1;
  }
  return result;
}

如您所见，它比前两个更快:

# cpp_vectorized ----------------------------------------------------------

cpp_vect <- function(bigg) {
  res_cppvect <- cppnextup(bigg, 0.5)
}

    # benchmark ---------------------------------------------------------------
    
    microbenchmark::microbenchmark(base_(),
                                   cpp_(), 
                                   cpp_vect(),
                                   times=3)
           expr       min        lq      mean    median        uq       max neval
        base_() 2014.7211 2016.8679 2068.9869 2019.0146 2096.1198 2173.2250     3
         cpp_()  663.0874  666.1540  718.5863  669.2207  746.3357  823.4507     3
     cpp_vect()  214.1745  221.2103  223.9532  228.2460  228.8426  229.4392     3

但当我在参数中传递一个更大的向量时，它就冻结了，永远不会返回结果.

res <- cpp_vect(bigg=rnorm(1000000)) # freezes

欢迎任何帮助.