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

Question

我在 R 中有一个循环，它非常慢（但有效）。目前，这个计算在我的笔记本电脑上大约需要 3 分钟，我认为它可以改进。最后，我将根据这段代码的结果循环遍历许多数据文件运行ning 计算，如果可能的话，我想使当前代码更快。

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

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

我正在寻找有关如何改进 R 代码的反馈and/or 更正 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) error causing an undefined behavior (UB)！您可以在将来通过将向量访问器从 [] 更改为 () 并将矩阵访问器从 () 更改为 .at() 来检测到这一点。

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

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

表示索引超出范围，因为索引必须始终小于范围（例如向量的长度 - 1）在 0 到 1 之间。

初步看来，这个问题主要是由于没有正确地将基于一的索引映射到基于零的索引引起的。

在使用 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 等价物的输出。

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

Optimizing loop on large data file in R, perhaps using Rcpp

performance

loops

r

rcpp