r - 优化 R 中大数据文件的循环，可能使用 Rcpp

Question

我在 R 中有一个循环，它很慢(但有效)。目前，这个计算在我的笔记本电脑上大约需要 3 分钟，我认为它可以改进。最终，我将遍历许多基于此代码结果运行计算的数据文件，如果可能，我希望使当前代码更快。

基本上，对于每个日期，对于 11 个不同的 X 值，循环获取过去 X 年的降雨值 (Y)，找到线性逆加权 (Z)，以便最旧的降雨值加权最小，乘以降雨(Y) 和权重 (Z) 得到一个向量 A，然后将 A 的和作为最终结果。这是针对数千个日期完成的。

但是，我想不出或找到任何建议来使 R 中的速度更快，因此我尝试在 Rcpp 中重写它，我对此知之甚少。我的 Rcpp 代码没有完全复制 R 代码，因为结果矩阵与它应该是的不同(错误)(out1 vs out2；我知道 out1 是正确的)。 Rcpp 代码似乎更快，但我只能使用几列对其进行测试，因为如果我尝试运行所有 11 列(i <= 10)，它就会开始崩溃(RStudio 中的 fatal error )。

我正在寻找有关如何改进 R 代码和/或更正 Rcpp 代码以提供正确结果而不是在过程中崩溃的反馈。

(虽然我在下面发布的代码没有显示它，但数据以 [作为数据框] 的方式加载到 R 中，用于在所示代码之外完成的一些计算。对于此处显示的特定计算，只有列使用了数据帧的 2。)

数据文件在这里:https://drive.google.com/file/d/0Bw_Ca37oxVmJekFBR2t4eDdKeGM/view?usp=sharing

R中的尝试

library(readxl)

library(readxl)
library(Rcpp)
file = data.frame(read_excel("lake.xlsx", trim_ws=T)
col_types=c("date","numeric","numeric","date",rep("numeric",4),"text")))
file[,1] = as.Date(file[,1], "%Y/%m/%d", tz="UTC")
file[,4] = as.Date(file[,4], "%Y/%m/%d", tz="UTC")

rainSUM = function(df){
rainsum = data.frame("6m"=as.numeric(), "1yr"=as.numeric(), "2yr"=as.numeric(), "3yr"=as.numeric(), "4yr"=as.numeric(), "5yr"=as.numeric(), "6yr"=as.numeric(), "7yr"=as.numeric(), "8yr"=as.numeric(), "9yr"=as.numeric(), "10yr"=as.numeric()) # create dataframe for storing the sum of weighted last d values

  Tdays <- length(df[,1])

  for(i in 1:11) {           # loop through the lags
    if (i==1) {
      d <- 183               # 6 month lag only has 183 days,
    } else {
      d <- (i-1)*366         # the rest have 366 days times the number of years
    }
    w <- 0:(d-1)/d

    for(k in 1:Tdays) {      # loop through rows of rain dataframe (k = row)
      if(d>k){               # get number of rain values needed for the lag
        rainsum[k,i] <- sum(df[1:k,2] * w[(d-k+1):d])                                
      } else{
        rainsum[k,i] <- sum(df[(k-d+1):k,2] * w)                          
      }
    }
  }
  return(rainsum)
}

out1 <- rainSUM(file)

尝试 Rcpp

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]

NumericVector myseq(int first, int last) {  // simulate R's X:Y sequence (step of 1)
  NumericVector y(0);
  for (int i=first; i<=last; ++i)
    y.push_back(i);
  return(y);
}

// [[Rcpp::export]]   
NumericVector splicer(NumericVector vec, int first, int last) {  // splicer
  NumericVector y(0);
  for (int i=first; i<=last; ++i)
    y.push_back(vec[i]);
  return(y);
}

// [[Rcpp::export]]              
NumericVector weighty(int d) {             // calculate inverse linear weight according to the number of days in lag 
  NumericVector a = myseq(1,d);            // sequence 1:d; length d
  NumericVector b = (a-1)/a;               // inverse linear
  return(b);                               // return vector                              
} 

// [[Rcpp::export]]              
NumericMatrix rainsumCPP(DataFrame df, int raincol) {
  NumericVector q(0);
  NumericMatrix rainsum(df.nrows(), 11);       // matrix with number of row days as data file and 11 columns 
  NumericVector p = df( raincol-1 );           // grab rain values (remember C++ first index is 0)
  for(int i = 0; i <= 10; i++) {                // loop through 11 columns (C++ index starts at 0!)
    if (i==0) {
      int d = 183;                             // 366*years lag days 
      NumericVector w = weighty(d);            // get weights for this lag series
      for(int k = 0; k < df.nrows(); k++) {    // loop through days (rows)
        if(d>k){                               // if not enough lag days for row, use what's available
          NumericVector m = splicer(p, 0, k);  // subset rain values according to the day being considered
          NumericVector u = splicer(w, (d-k), (d-1));   // same for weight
          m = m*u;                              // multiply rain values by weights
          rainsum(k,i) = sum(m);               // add the sum of the weighted rain to the rainsum matrix
        } else{
          NumericVector m = splicer(p, k-d+1, k);
          m = m*w;
          rainsum(k,i) = sum(m);
        }
      }
    }
    else {
      int d = i*366;                           // 183 lag days if column 0
      NumericVector w = weighty(d);            // get weights for this lag series
      for(int k = 0; k < df.nrows(); k++) {    // loop through days (rows)
        if(d>k){                               // if not enough lag days for row, use what's available
          NumericVector m = splicer(p, 0, k);  // subset rain values according to the day being considered
          NumericVector u = splicer(w, (d-k), (d-1));   // same for weight
          m = m*u;                             // multiply rain values by weights
          rainsum(k,i) = sum(m);               // add the sum of the weighted rain to the rainsum matrix
        } else{
          NumericVector m = splicer(p, k-d+1, k);
          m = m*w;
          rainsum(k,i) = sum(m);
        }
      } 
    }
  }
  return(rainsum);
}


/*** R
out2 = rainsumCPP(file, raincol) # raincol currently = 2
  */

Answer 1

恭喜!您有一个 index out of bounds (OOB)导致 undefined behavior (UB) 的错误!您可以在将来通过从 [] 更改矢量访问器来检测到这一点。至 ()以及来自 () 的矩阵访问器(accessor)至 .at() .

切换到这些访问器会产生:

Error in rainsumCPP(file, 2) : 
  Index out of bounds: [index=182; extent=182].

这表明索引越界，因为索引必须始终小于范围(例如向量的长度 - 1)在 0 到 1 之间。

初步观察表明，此问题主要是由于未将基于 1 的索引正确映射到基于 0 的索引造成的。

在玩弄 myseq() 时, splicer() , 和 weighty()函数，它们与给定的 R 等效输入不匹配。这可以通过使用 all.equal(R_result, Rcpp_Result) 来检查。 .这种不匹配分为两部分:1. 两者的界限 myseq和 splicer和 2. 内部完成的反转 weighty .

因此，通过使用以下修改后的函数，您应该为获得正确的结果打下良好的基础。

// [[Rcpp::export]]
NumericVector myseq(int first, int last) {  // simulate R's X:Y sequence (step of 1)
  int vec_len = abs(last - first);

  NumericVector y = no_init(vec_len);

  int count = 0;
  for (int i = first; i < last; ++i) {
    y(count) = count;
    count++;
  }

  return y;
}

// [[Rcpp::export]]   
NumericVector splicer(NumericVector vec, int first, int last) {  // splicer

  int vec_len = abs(last - first);

  NumericVector y = no_init(vec_len);

  int count = 0;
  for (int i = first; i < last; ++i) {
    y(count) = vec(i);

    count++;
  }

  return y;
}

// [[Rcpp::export]]              
NumericVector weighty(int d) {       // calculate inverse linear weight according to the number of days in lag 
  NumericVector a = myseq(0, d - 1); // (fixed) sequence 1:d; length d
  NumericVector b = a / d;           // (fixed) inverse linear
  return(b);                         // return vector                              
}

从那里，您可能需要修改 rainsumCpp因为没有给出 R 等价物的输出。