おそらくRcppを使って、Rの大きなデータファイルのループを最適化します。

Question

私は非常に遅いですが（しかし動作する）Rのループを持っています。現在のところ、この計算にはノートパソコンで約3分かかりますが、改善できると思います。最終的には、このコードの結果に基づいて計算を実行している多くのデータファイルをループします。できるだけ早く現在のコードを作成したいと思います。

基本的に、各日付について、11個の異なる値のXについて、ループは最後のX年分の降雨値（Y）を取得し、最も古い降雨値が最小になるように線形逆重み付け（Z）ベクトルAを得るために雨（Y）と重み（Z）を掛け合わせ、最後の結果としてAの合計をとる。これは何千もの日付で行われます。

しかし、Rでこれをもっと速くする方法はないと思うので、Rcppに書き直そうとしました。そこでは私は知識が限られています。私のRcppコードは、Rコードを正確に複製しません。なぜなら、結果として得られる行列は、（out1対out2;私はout1が正しいことを知っています。 Rcppのコードは高速だと思われますが、11列すべてを実行しようとすると（RSTudioの致命的なエラー）、それがクラッシュし始めるため、いくつかの列を使ってしかテストできません（i < = 10）。

私はRコードを改善したり、正しい結果を提供し、プロセスでクラッシュしないようにRcppコードを修正する方法に関するフィードバックを探しています。

（以下に掲載したコードでは表示していませんが、表示されているコードの外で行われたいくつかの計算では、データフレームとしてデータが読み込まれます。ここで、データフレームの唯一のカラム2が使用される）

データファイルはここにある：。RcppにおけるRでhttps://drive.google.com/file/d/0Bw_Ca37oxVmJekFBR2t4eDdKeGM/view?usp=sharing

試み

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) 

試み

#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)が発生しています。ベクトルアクセサーを[]から()に変更し、マトリックスアクセサーを()から.at()に変更することで、これを検出できます。

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

インデックスは常にエクステント（ - 1ベクターの例えば長さ）よりも0及び1未満の間でなければならないようにインデックスが範囲外であることを示している：これらのアクセサ収量への切り替え

。

この問題は、1ベースのインデックスをゼロベースのインデックスに正しくマッピングしないことが主な原因であることを示しています。 myseq()、splicer()、およびweighty()機能で遊んでたら

、彼らはは彼らのR同等与えられた入力と一致していません。これはall.equal(R_result, Rcpp_Result)を使用して確認できます。このミスマッチは2つの部分に分かれています：1. myseqとsplicerの両方の境界と、内部で行われた反転が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        
} 

何も出力がR同等だったもののために与えられなかったとして、あなたはおそらくrainsumCppを変更する必要があります。