Rcppで不等間隔の時系列データに対して移動平均を計算する
一年ぶりの更新ですね…。
不等間隔の時系列データに対して一定のwindow幅で移動平均や移動分散などを計算したいが、Rでは時間がかかる。という相談が slack の r-wakalang にありました。
推測ですが、Rで普通に（?）やるとwindowの範囲に入っている値を検索するのに、毎回全てのデータから検索してしまうところがボトルネックになっていると思います。そこでRcppを使って値の検索の部分を効率化してみました。
今回は移動平均を計算していますが、紹介する関数を雛形にすれば、windowに対して様々な処理を適用できるはずです。
まずはデータを準備します。
# windowの始点 from <- seq.POSIXt(as.POSIXlt("2017-01-01 00:00:00"), as.POSIXlt("2017-01-31 00:00:00"), by = 60*60) # windowの終点 to <- from + 3*60*60 # 不等間隔な計測データを適当に再現する library(dplyr) time <- seq.POSIXt(min(from), max(to), by = 10) %>% sample(size = 10000) %>% sort input <- data.frame( time, #計測時刻 data = runif(length(time), 0, 5) #計測データ )
こんな感じのデータ。
> input
time data
1 2017-01-01 00:00:00 4.573989e+00
2 2017-01-01 00:02:40 1.899510e+00
3 2017-01-01 00:06:20 3.148162e+00
4 2017-01-01 00:10:20 1.928675e+00
5 2017-01-01 00:11:30 2.711351e+00
6 2017-01-01 00:13:30 2.464259e+00
7 2017-01-01 00:15:50 7.405586e-01
8 2017-01-01 00:16:40 2.433681e+00
9 2017-01-01 00:19:30 4.152487e+00
...略
Rバージョン
まずは、Rバージョン。工夫すればもっと効率化できるかもしれないですが、とりあえず最初の思いつきで書いたやつ。
# Rで移動平均（デフォルトのwindow幅は3時間） moving_mean_r <- function(input, from, time_window=3*60*60){ to <- from + time_window mean <- mapply(function(f, t){input %>% filter(time>=f, time<t) %>% '[['("data") %>% mean},from,to) data.frame(from, to, mean) }
上のコードがあまりにも遅いので書き直しました。（ここでdplyr::filter使ってはいけなかったんや…。）
moving_mean_r2 <- function(input, from, time_window=3*60*60){ to <- from + time_window mean <- mapply(function(f, t){input[input$time>=f & input$time<t,"data"] %>% mean;} ,from,to) data.frame(from, to, mean) }
Rcppバージョン
#include <Rcpp.h> #include <algorithm> using namespace Rcpp; //[[Rcpp::plugins(cpp11)]] // [[Rcpp::export]] DataFrame moving_mean_rcpp( DataFrame input, //（計測時刻time, 計測値data）のデータフレーム newDatetimeVector from, //移動平均を計算するタイミングの始点ベクトル double time_window = 3*60*60) //移動平均を計算するwindow幅（秒） { // 計測時刻と計測値をベクトルとして取り出す newDatetimeVector time = input["time"]; // 今回は time は昇順にソートされているのが前提です。 NumericVector data = input["data"]; // 移動平均を計算するタイミングの終点ベクトル newDatetimeVector to = from + time_window; // 移動平均を計算する数 R_xlen_t N = from.length(); // 移動平均を格納するベクトル NumericVector value(N); // ベクトル要素の位置をあらわすオブジェクト newDatetimeVector::iterator begin = time.begin(); newDatetimeVector::iterator end = time.end(); newDatetimeVector::iterator p1 = begin; newDatetimeVector::iterator p2 = begin; // window i についてループ for(R_xlen_t i = 0; i < N; ++i){ // Rcout << "i=" << i << "\n"; double f = from[i]; //windowの始点の時刻 double t = f + time_window; //windowの終点の時刻 // windowの終点が最初の計測時刻以前の時はNA、または // windowの始点が最後の計測時刻のより後の時はNA if(t <= *begin || f > *(end-1)){ value[i] = NA_REAL; continue;//次のループへ } // ベクトル time の位置 p1 以降の要素xから // 時刻がwindowの始点f「以降」である「最初の要素」の位置を p1 とする p1 = std::find_if(p1, end, [&f](double x){return f<=x;}); // p1 = std::lower_bound(p1, end, f); //上と同義 // ベクトル time の位置 p1 以降の要素xから // 時刻がwindowの終点t「より前」である「最後の要素」の位置を p2 とする // （下では、時刻がwindowの終点t「以降」である「最初の要素」の１つ前の位置、にすることで実現している’） p2 = std::find_if(p1, end, [&t](double x){return t<=x;}) - 1 ; // p2 = std::lower_bound(p1, end, t) - 1 ;//上と同義 // 要素の位置p1,p2を、要素番号i1, i2に変換する R_xlen_t i1 = p1 - begin; R_xlen_t i2 = p2 - begin; // 要素番号の確認 // C++は要素番号が0から始まるのでRに合わせるために1を足している // Rcout << "i1 = " << i1+1 << " i2 = " << i2+1 << "\n"; // 該当する範囲のデータについて平均を計算する if(i1>i2) { value[i] = NA_REAL; // window内にデータがない場合 } else { value[i] = mean(data[seq(i1, i2)]); } // ↑を変更することで様々なwindow関数を作成できる } // 平均を計算した時間の範囲と、平均の値をデータフレームとして出力する DataFrame out = DataFrame::create( Named("from", from), Named("to", to), Named("mean", value)); return out; }
RバージョンとRcppバージョンの結果の一致を確認
> # Rバージョン > res_r <- moving_mean_r(input, from) > > # Rcppバージョン > res_rcpp <- moving_mean_rcpp(input, from) > > # 結果の比較（数値計算の誤差も考慮して） > isTRUE(all.equal(res_r$mean, res_rcpp$mean)) [1] TRUE
パフォーマンス比較
> microbenchmark::microbenchmark(
+ moving_mean_r(input, from),
+ moving_mean_r2(input, from),
+ moving_mean_rcpp(input, from)
+ )
Unit: microseconds
expr min lq mean median uq max neval
moving_mean_r(input, from) 904298.672 951065.9675 985553.3977 961292.6805 1015783.925 1179497.694 100
moving_mean_r2(input, from) 191614.227 208638.5365 226343.9220 216555.8075 230387.330 317652.639 100
moving_mean_rcpp(input, from) 441.276 488.5045 569.5601 510.9075 546.498 2998.015 100
RとRcppでかなり差が出ました。
Rcpp で疎行列をデータフレームに変換する
しばらく体調を崩してしまっていたため、すっかりご無沙してしまいました。ようやく回復してきましたので活動再開しようと思います。
というわけで、だいぶ間があいてしまいましたが、前回の続きということで、Matrixパッケージで作成した疎行列をデータフレームに変換する処理を紹介します。（正直、必要性はあまりないとは思いますが、Rcpp の使い方を示す例としては悪くないと思います。あと、疎行列の実装方法の一つである圧縮列格納方式については、説明する前に力尽きました。）
ちなみに、R の場合には以下のようにします。
# 疎行列の作成
library(Matrix)
sm <- Matrix(sample(c(0,1), size = 25, replace = TRUE)
, nrow = 5
, sparse = TRUE)
# データフレームへの変換
df <- as.data.frame(as.matrix(sm))
これを Rcpp で実装することで、以下の Rcpp 機能の使用例を示します。
- S4 クラスのスロットへのアクセス
- stop() によるエラー処理
- 要素数や要素番号の型としての R_len_t, R_xlen_t
- diff(), seq(), rep() など R と同じ機能の関数の利用
- NULL値の判定
- String 型文字列の結合（+=）
- サイズを指定したデータフレームの作成
- リストからデータフレームへの変換
（追記：2016/06/19）64 bit システムでは R のベクトルの要素数の最大値は理論上 「2 ^ 64 - 1」 なので、int の最大値「2 ^ 31 - 1」よりも大きくなり得ます。なので、基本的には要素数や要素番号の型としては R_len_t（int型と同じ）ではなく R_xlen_t （ptrdiff_t 型と同じ）を使うべきです。しかし、今回のケースで言うと sm.slot("Dim") は IntegerVector （要素の型は int ）になっているので、変数 nrow, ncol などは R_len_t でも良いです。
実装例
// [[Rcpp::export]] List sparseMatrixToDataFrame(S4 sm){ // Matrixパッケージの疎行列（dgCMatrix）ではないデータが // 渡された場合は拒否します if(as<std::string>(sm.attr("class"))!="dgCMatrix") stop("Input must be dgCMatrix."); // 行数と列数 //要素数や要素番号の型には R_xlen_t を使います IntegerVector dim = sm.slot("Dim"); R_xlen_t nrow = dim[0]; R_xlen_t ncol = dim[1]; // dgCMatrix から値を取り出します // dgCMatrix の疎行列の形式は圧縮列格納方式となっています NumericVector X = sm.slot("x");//非ゼロ要素の値 IntegerVector I = sm.slot("i");//非ゼロ要素の行番号（0ベース） IntegerVector P = sm.slot("p");//非ゼロ要素の列番号の開始位置（0ベース） // P から非ゼロ要素の列番号 J (0ベース) を再生します IntegerVector J(X.length()); IntegerVector repeat = diff(P); //各列番号が連続する数 R_xlen_t n = repeat.length(); for(R_len_t j=0; j<n; ++j){ if(repeat[j]==0) continue; //非ゼロ値が存在しない列は飛ばします // J に 各非ゼロ値に対応する列番号 j をセットします J[seq(R_xlen_t(P[j]),R_xlen_t(P[j+1]-1))] = rep(j, repeat[j]); } // 出力用データを作成します // ここでは列数が ncol のデータフレームを作成したいですが // DataFrame df(ncol)；とすると、１列目の値が ncol である // データフレームができてしまいます、。 // そこで要素数が ncol であるリストを作成してから // 最後にデータフレームに変換します List L(ncol); // リストを初期化するため // 要素の値が 0 で長さが nrow のベクトルを作成して // そのベクトルをリストの各要素として代入します for(R_xlen_t i=0;i<ncol;++i){ L[i] = NumericVector(nrow); } // 出力用データに値をセットします // 行番号 I、列番号 J、の要素に値 X を代入します n = I.length(); NumericVector column; for(R_xlen_t k=0; k<n; ++k){ // column は L の J[k] 番目の要素（列）への参照となります column = L[J[k]]; // column の I[k] 番目の要素に X[k] の値を代入します column[I[k]] = X[k]; } // 出力用データにセットする行名・列名を作成します List Dimnames = sm.slot("Dimnames"); CharacterVector rownames(nrow); CharacterVector colnames(ncol); // 行名 if(Dimnames[0]==R_NilValue){ //行名がNULLの場合は、1 2 ... という行名にします IntegerVector v = seq_len(nrow); rownames = as<CharacterVector>(v); } else { rownames=Dimnames[0]; } // 列名 if(Dimnames[1]==R_NilValue){ //列名がNULLの場合は、V1 V2 ... という列名にします //R の paste() 関数を呼んだら遅かったのでやめました for(int i=0; i<ncol; ++i){ colnames[i] = (String("V") += String(i+1)); } } else { colnames=Dimnames[1]; } // これらの属性に値をセットして返すと // R ではデータフレームとして扱われます L.attr("row.names") = rownames; L.attr("names") = colnames; L.attr("class") = "data.frame"; return L; }
実行結果
> sm <- Matrix(sample(c(0,1), size = 25, replace = TRUE)
+ , nrow = 5
+ , sparse = TRUE)
> sm
5 x 5 sparse Matrix of class "dgCMatrix"
[1,] 1 1 1 1 .
[2,] . 1 . 1 .
[3,] . . 1 1 1
[4,] 1 1 1 1 .
[5,] . 1 1 1 1
> sparseMatrixToDataFrame(sm)
V1 V2 V3 V4 V5
1 1 1 1 1 0
2 0 1 0 1 0
3 0 0 1 1 1
4 1 1 1 1 0
5 0 1 1 1 1
パフォーマンス比較
それでは R バージョンと速度を比較してみましょう
> sm <- Matrix(sample(c(0,1), size = 1000000, replace = TRUE)
+ , ncol = 10000
+ , sparse = TRUE)
> microbenchmark::microbenchmark(sparseMatrixToDataFrame(sm)
+ ,as.data.frame(as.matrix(sm))
+ ,times = 10
+ )
Unit: milliseconds
expr min lq mean median uq max neval
sparseMatrixToDataFrame(sm) 38.63905 45.4774 77.11882 55.54285 62.58709 288.8416 10
as.data.frame(as.matrix(sm)) 129.57834 140.7170 258.36984 261.72998 374.73649 391.0923 10
中央値でみるとRcpp版は 55.5 ミリ秒に対して、R版は 261.7 ミリ秒なので、Rcpp のほうが 1/5 の時間で終わっています。実際のところR版も関数はC言語で実装されているので結構速くて抜くのに苦労したのですが、おかげで有用な知見が得られました。
Rcppでデータフレームを疎行列に変換する
前回までは、Rcpp を使ってファイルの行数をカウントするコードを紹介しました。しかし、紹介したコードでは C/C++ の関数やクラスしか使わなかったため、Rcpp の機能が活用されていない感じになってしまいました。そこで今回は Rcpp の機能を使って、いかに R との連携がスムーズに行われるか見て見たいと思います。その例として、データフレームを疎行列に変換するプログラムを紹介します。
今回紹介するプログラムは Rcpp がもつ R との連携機能を色々活用しています。
- DataFrame の要素へのアクセス
- DataFrame の行数や列数
- NA 値の判定
- パッケージの関数の呼び出し
- オブジェクトの属性値（列名、行名）へのアクセス
- NULL値のセット
- 標準 C++ のデータ構造を利用
コード例
まずは、さっそくコード例を見てみましょう。ちなみに、ここで紹介しているプログラムでは、カテゴリカル変数のダミー変数への変換は、既に行われているという前提になっておりますので、あしからず。
// [[Rcpp::export]] S4 asSparseMatrix( DataFrame df, bool na = 0){ // df : データフレーム、値は数値のみである前提 // na df にNAが含まれる場合は 1、 含まれていない場合は 0 にセット //データの行数と列数 int nrow = df.nrows(); int ncol = df.length(); // R の Vector では要素の追加（push_back）の効率が悪いので // 標準 C++ の vector を使用する。 // R_xlen_t は行番号や列番号の型 std::vector<R_xlen_t> rows; std::vector<R_xlen_t> cols; std::vector<double> vals; //データフレームの全ての要素にアクセスして、 //0じゃない要素の位置と値を保存する if(na){//NAありの場合 //データがNAを含む可能性がある場合は //NAのチェックを行う for(R_xlen_t col=0; col<ncol; ++col){ NumericVector column = df[col]; for(R_xlen_t row=0; row<nrow; ++row){ // NA の要素は 0 とみなす if(NumericVector::is_na(column[row])) continue; else if(column[row]!=0.0){ rows.push_back(row+1); cols.push_back(col+1); vals.push_back(column[row]); } } } } else { //データがNAを含まないと保証された場合は //NAのチェックを行わない for(R_xlen_t col=0; col<ncol; ++col){ NumericVector column = df[col]; for(R_xlen_t row=0; row<nrow; ++row){ if(column[row]!=0.0){ rows.push_back(row+1); cols.push_back(col+1); vals.push_back(column[row]); } } } } // Matrix パッケージから sparseMatrix 関数を読み込み // Matrix::sparseMatrix() という形式で呼び出したのと同じ Environment env = Environment::namespace_env("Matrix"); Function sparseMatrix = env["sparseMatrix"]; // wrap() を使って std::vector を NumericVector に変換している。 S4 sm = sparseMatrix( Named("i") = wrap(rows), Named("j") = wrap(cols), Named("x") = wrap(vals), Named("dims") = NumericVector::create(nrow,ncol)); //疎行列に列名をセットする、行名は NULL にする List dimnames = List::create(R_NilValue, df.names()); sm.attr("Dimnames") = dimnames; return sm; }
実行例
x <- sample(c(0,1), size = 5, replace = TRUE) y <- sample(c(0,1), size = 5, replace = TRUE) data <- data.frame(x,y) sm1 <- Matrix::Matrix(as.matrix(data), sparse = TRUE) sm2 <- asSparseMatrix(data) identical(sm1, sm2) // TRUE
この関数を作成した背景
Matrix パッケージの疎行列クラス（dgCMatrix）を分析で使用しているのですが、データフレームから疎行列を作成するためには、通常はMatrix(as.matrix(data), sparse = TRUE) という感じで疎行列に変換します。しかし、データが巨大になってくると、行列のサイズの限界により as.matrix() での行列変換ができなくなってしまうのです。それならデータフレームから疎行列を直接変換すればいいじゃん、というのが自然だと思うのですが、残念ながら Matrix パッケージでは対応してくれておりません。
実は、Matrixパッケージには、疎行列を作成するもう一つの関数 sparseMatrix() があり、こちらは全ての非ゼロの要素の行番号・列番号・値を入力することで疎行列を作成することができます。そのためにはデータの全ての要素にアクセスして値が、ゼロかどうか判定し、行番号・列番号と値を記録しなければなりませんが、このように値をいちいちチェックする処理は R の苦手とする処理の１つです。そこで Rcpp の出番というわけです。
標準 C++ データ構造を使用する理由
コード例では、非ゼロの要素の行番号・列番号・値を記録するために、NumericVector の代わりに、標準 C++ のデータ構造 std::vector を用いています。データフレームの中に非ゼロ要素がいくつあるかは、全ての要素をチェックしてからでないとわからないので、あらかじめ決められたサイズのベクトルに値を保持することが難しいです。そこで、ベクトルのサイズを動的に増加させることができると嬉しいのですが、Rcpp の Vector（というか R のベクトル）は、一旦ベクトルを作成した後で、要素のサイズを変更する （v <- c(v, 1)）処理を行うと、ベクトルの要素の値を全てをコピーする処理が発生するので、非常に効率が悪いという問題があります。しかし、標準 C++ の std::vector は、新たな要素の追加の効率に優れていますので、こちらを使用しました。
パフォーマンス
Rcpp のコード例のプログラムのパフォーマンスを評価するため、比較対象として、Matrix(as.matrix(data), sparse = TRUE) の他に、コード例と同じ処理を R で実装したバージョンを用意しました。今時なら dplyr パッケージを利用すればもっと効率がよい実装が可能ですが、Rの基本機能で実装すると以下のようにになると思います。
asSparseMatrix_R <- function(data){
tmp <-
lapply(1:ncol(data), function(col){
row <- which(data[,col]!=0)
data.frame(row=row, col=col, val=data[row,col])
})
res <- do.call("rbind", tmp)
Matrix::sparseMatrix(i = res[,1],
j = res[,2],
x = res[,3],
dims = c(nrow(data),ncol(data)))
}
それではパフォーマンスを比較してみましょう。
> nrow <- 10000
> ncol <- 100
> m <- matrix(sample(c(0,1), size = nrow*ncol, replace = TRUE)
+ , nrow = nrow)
> data <- as.data.frame(m)
>
> microbenchmark::microbenchmark(
+ Matrix::Matrix(as.matrix(data), sparse = TRUE),
+ asSparseMatrix(data, na = 0),
+ asSparseMatrix_R(data),
+ times = 10)
Unit: milliseconds
expr min lq mean median uq max neval
Matrix::Matrix(as.matrix(data), sparse = TRUE) 35.88858 41.0130 97.69263 49.57377 58.56511 306.7089 10
asSparseMatrix(data) 113.41807 127.1215 159.82937 134.94027 150.90175 391.1737 10
asSparseMatrix_R(data) 863.00960 1104.5381 1248.25324 1184.21548 1485.64696 1579.7606 10
残念ながら、この比較では Matrix() 関数に負けてしまいましたが、まずまずの結果です。もともとMatrix() 関数が使えない状況のためにasSparseMatrix() は作成したので、これで十分でしょう。ここで使用したデータはサイズが7MBくらいのデータですが、数十GB程度のデータに対して、それほど時間もかからず処理できるので、業務でも役に立ってくれています。
というわけで、今回は以上です。いかがでしたでしょうか、ちょっと今回は機能が盛り沢山な例になってしまったので、消化するのが大変かもしれないですね。Rcppの様々な機能の使い方が端的に示されていますので、ゆるゆると学んでいただけば良いかなと思います。
次回は、今回とは逆に疎行列をデータフレームに変換する処理を紹介する予定です。
Rcpp でファイルの行数をカウントする（２）
前回に引き続きファイルの行数のカウントに取り組みます。
前回は getline() と istreambuf_iterator() を使った方法を紹介しましたが、残念ながら Unix コマンド wc には負けてしまいました。これをこのままにしておくのも口惜しいので、C言語で書かれた wc のソースコードからファイル行数をカウントする部分を参考にして、それを関数として実装しなおしました。
ソースを見ると wc コマンドではファイルを決められた文字数ずつメモリに読み込みながら行数をカウントしていました。そこで、この一度に読み込む文字数をユーザーが変えられるように Rcpp の関数にしてみました。
rcpp_wc.cpp
#include <Rcpp.h> using namespace Rcpp; // [[Rcpp::export]] unsigned int rcpp_wc(std::string file, int buf_size=10000) { //file : ファイルへのパス //buf_size : 一度の読みこむ文字数 //ファイルを開く FILE *fp; fp = fopen (file.c_str(), "r"); uintmax_t lines=0;//行数のカウント用 char buf[buf_size];//読み込んだ文字を保存する変数 size_t bytes_read;//実際に読み込まれたバイト数（文字数）を格納する変数 //行が長い場合には、行の数え方を変えるためのフラグ bool long_lines = false; //ファイルをから buf_size ずつ文字列を読み込んでゆく while ((bytes_read = fread(buf, 1, buf_size, fp)) > 0) { char *p = buf; //ポインタ p にバッファーの先頭文字のアドレスを与える char *end = p + bytes_read; //このステップで読み込まれた文字列の最後の文字へのアドレスを保持する uintmax_t plines = lines; //前のステップで読み込んだ文字列中にあった改行の数を plines に保存しておく // 改行 '\n' の数を数える if (! long_lines)//1行が短い場合 {//文字列を1文字ずつ走査して改行をカウント while (p != end) lines += *p++ == '\n'; } else {//1行が長い場合 // memchr()の呼び出しにかかる時間を考慮しても memchr() を使ったほうが速度が速い // memchr(p, '\n', n) は位置 p にある文字から n 文字先のまでの範囲から '\n' を探して // 最初に見つけた '\n' の位置を返す関数 // reinterpret_cast<char*>() は memchr() の返値の型である void* を char* として解釈するように命令している while ((p = reinterpret_cast<char*>(memchr (p, '\n', end - p)))) // 最初に見つかる '\n' を探す { ++p;//見つけた '\n' の次の文字に移る ++lines;//見つけた '\n' の数（行数）をカウントアップ } } //バッファー中に含まれる行数の平均値が 15 以上になったら、 //次のステップからは memchr() を使うようにフラグを立てる if (lines - plines <= bytes_read / 15) long_lines = true; else long_lines = false; //次のステップへ } return lines; }
では、これを元の wc コマンドと比較してみましょう。一度に読み込む文字数を変えると行数のカウントのスピードはどう変わるでしょうか。
sourceCpp("rcpp_wc.cpp")
microbenchmark::microbenchmark(
system("wc -l hoge.csv"),
rcpp_wc("hoge.csv",10000),
rcpp_wc("hoge.csv",1000),
rcpp_wc("hoge.csv",100),
rcpp_wc("hoge.csv",10),
times=5)
実行結果
Unit: milliseconds
expr min lq mean median uq max neval
system("wc -l hoge.csv") 1340.7916 1403.8807 1735.8336 1498.6164 2183.7768 2252.1026 5
rcpp_wc("hoge.csv", 10000) 345.6389 382.7236 476.2079 388.7028 423.4235 840.5509 5
rcpp_wc("hoge.csv", 1000) 601.1652 618.4181 675.9131 632.3813 655.2329 872.3680 5
rcpp_wc("hoge.csv", 100) 1725.2825 1728.1933 1799.3132 1816.1720 1820.8418 1906.0762 5
rcpp_wc("hoge.csv", 10) 11762.0427 11845.1454 12273.0468 12159.5957 12731.1063 12867.3437 5
おおっ！ついに wc に勝利しました。一度に読み込む文字数を増やすと行数カウントのスピードも伸びていることがわかります。かなり劇的な効果です。hoge.csv は600MB程度のサイズですが、この程度なら瞬殺できることがわかります。こんなに速くなるとは思いませんでした。テストマシンがしょぼいのでこれ以上大きなファイルでは試してないですが、かなり大きなファイルでもこれは期待できそうです。
いやー、これでまた業務が捗ってしまいますね（ドヤ）。
ということで、Rcpp の凄さの一端は感じていただけたかと思いますが、コードが Rcpp というより C/C++ そのまんまな感じになってしまったので、次回はもうちょっと Rcpp の機能を利用したコードを紹介したいと思います。
でわー
Rcpp でファイルの行数をカウントする
新しいサイトでブログを再開してみました。名前はそのうち考えます。
とりあえず、せっかく Rcpp 入門を GitBook でまとめたところなので、布教のために Rcpp の活用例をちまちまと紹介して行こうかなと思います。
記念すべき最初の記事として teramonagi さんの 「ファイルの行数を調べたい」 に応えてみました。
ここでは C++ を用いた２通りの実装を紹介します。 よくみたら Rcpp のクラス全然使ってないですね。でも、これで OK なのです。
Rcpp コード
count_lines.cpp
#include <Rcpp.h> #include <fstream> // [[Rcpp::export]] unsigned int Rcpp_count_lines1(std::string path){ std::ifstream inFile(path.c_str()); return std::count(std::istreambuf_iterator<char>(inFile), std::istreambuf_iterator<char>(), '\n'); } // [[Rcpp::export]] unsigned int Rcpp_count_lines2(std::string path){ unsigned int n = 0; std::string line; std::ifstream myfile(path); while (std::getline(myfile, line)) ++n; return n; }
では、さっそく teramonagi さんの記事を参考に wc、 R.utils::countLines( ) と比較してみましょう！
Rcpp::sourceCpp('count_lines.cpp')
microbenchmark::microbenchmark(
system("wc -l hoge.csv"),
R.utils::countLines("hoge.csv"),
Rcpp_count_lines1("hoge.csv"),
Rcpp_count_lines2("hoge.csv"),
times=5)
実行結果
Unit: seconds
expr min lq mean median uq max neval
system("wc -l hoge.csv") 1.225005 1.240038 1.519176 1.317884 1.436273 2.376679 5
R.utils::countLines("hoge.csv") 13.437869 14.185085 15.571043 15.480933 16.912233 17.839095 5
Rcpp_count_lines1("hoge.csv") 2.959419 3.066835 3.308570 3.261107 3.618280 3.637206 5
Rcpp_count_lines2("hoge.csv") 8.088860 8.184448 8.447557 8.358582 8.715715 8.890180 5
うーん、残念！ wc には負けてしまいました。さすがに wc 速いですね。でも R.utils::countLines() には余裕で勝ちました。C++ でも std::getline() より std::istreambuf_iterator を使った方法のほうが倍以上高速なようです。
ただ wc と Rcpp_count_lines1 は改行コードを数えているだけなので、最後の行のラストに改行がない場合には、行数が１つ少なくカウントされてしまうのが、ちと残念ですね。Rcpp_count_lines2 の getline の方はちゃんとカウントしてくれます。
C++ というと機械学習アルゴリズムの実装のような高度な用途にだけ使うものだろうと思われるかもしれないですが、今回の例のように、ちょっとした日常業務のツールを作成する場面でも案外役に立ってくれそうです。
というわけで最初の記事は以上です。いかがでしたでしょうか？ これからこのような感じで Rcpp の活用例を紹介してゆきたいと思います。
ではでは〜。
