[[Rの備忘録]] [[dirk.eddelbuettel Rcpp:http://dirk.eddelbuettel.com/code/rcpp.html]] #contents * 基本 [#hbf6ee84] Rcppインターフェイスを使ったクラスライブラリを作成するコード側で,たとえば Rから受け取った SEXP 型 a は,以下のように Rcpp::RObjectクラスの派生クラスに変換する必要があるのだが, Rcpp::NumericVector xa(a); この xa はそのまま演算できる. #code(c){{ #highlighter(CPP,number){{ RcppExport SEXP convolve3cpp(SEXP a, SEXP b) { Rcpp::NumericVector xa(a); Rcpp::NumericVector xb(b); int n_xa = xa.size(), n_xb = xb.size(); int nab = n_xa + n_xb - 1; Rcpp::NumericVector xab(nab); for (int i = 0; i < n_xa; i++) for (int j = 0; j < n_xb; j++) xab[i + j] += xa[i] * xb[j]; return xab; } }} これをRに返す際には,明示的に型変換する必要がないので,ついつい SEXP そのものだと誤解しそうになるが,あくまで Rcpp::RObjectクラスのオブジェクトであることを忘れないようにしないと,妙な代入などをおこなってしまいそうになる. ちなみに,伝統的な R API で,Rinternals.h を使う場合は,演算処理のためにポインタをとる. #code(x){{ #highlighter(CPP,number){{ SEXP convolve2(SEXP a, SEXP b) { R_len_t i, j, na, nb, nab; double *xa, *xb, *xab; SEXP ab; // 型を指定して保護 PROTECT(a = coerceVector(a, REALSXP)); PROTECT(b = coerceVector(b, REALSXP)); na = length(a); nb = length(b); nab = na + nb - 1; PROTECT(ab = allocVector(REALSXP, nab)); // C言語での演算のためポインタを取得 xa = REAL(a); xb = REAL(b); xab = REAL(ab); for(i = 0; i < nab; i++) xab[i] = 0.0; for(i = 0; i < na; i++) for(j = 0; j < nb; j++) xab[i + j] += xa[i] * xb[j]; UNPROTECT(3); //SEXP型なのでそのまま返却できる return(ab); } }} Rdefines.h を使う場合も,演算処理のためにポインタをとる. #code(c){{ #highlighter(php,number){{ SEXP convolve2(SEXP a, SEXP b) { R_len_t i, j, na, nb, nab; double *xa, *xb, *xab; SEXP ab; // 数値型として保護し PROTECT(a = AS_NUMERIC(a)); PROTECT(b = AS_NUMERIC(b)); na = LENGTH(a); nb = LENGTH(b); nab = na + nb - 1; PROTECT(ab = NEW_NUMERIC(nab)); // Cの演算のためにポインタを取り出す xa = NUMERIC_POINTER(a); xb = NUMERIC_POINTER(b); xab = NUMERIC_POINTER(ab); // ポインタを通した計算 for(i = 0; i < nab; i++) xab[i] = 0.0; for(i = 0; i < na; i++) for(j = 0; j < nb; j++) xab[i + j] += xa[i] * xb[j]; UNPROTECT(3); // SEXP型なのでそのまま返却できる return(ab); } }} と明示的にPROTECTが必要になる. * 引数を受ける [#a24d18a1] こんな感じで利用するのが一番楽? SEXP anyFunc ( SEXP arg1, SEXP arg2) { int a1 = as<int> (arg1); int a2 = as<int> (arg2); ... } 推奨は Rcpp::List rparam(params); std::string str = Rcpp::as<std::string>(rparam["moji"]); double num1 = Rcpp::as<double>(rparam["num1"]); int int1 = Rcpp::as<int>(rparam["num2"]); params <- list (moji ='hogehoge', num1= 3.14, num2 =10) func <- function ( parL ) { .Call("RcppParams", params = parL )} x <- func (params) でちとメンドイ 古い Rcpp (RcppClassic) では RcppParams par(params); double num1 = par.getDoubleValue("num1"); int int1 = par.getIntValue("num2"); string str = par.getStringValue("txt"); bool yn = par.getBoolValue("yesno"); * DataFrame 作成 [#x092ad8e] #include "Rcpp.h" RcppExport SEXP testDF (){ Rcpp::IntegerVector v = Rcpp::IntegerVector::create(1,2,3); std::vector<std::string> str(3); str[0] = "A"; str[1] = "B"; str[2] = "C"; return Rcpp::DataFrame::create(Rcpp::Named("Var1", v), Rcpp::Named ("Var2", s), Rcpp::Named("stringsAsFactors", false )); // あるいは return Rcpp::DataFrame::create(Rcpp::Named("VAR1") = v, Rcpp::Named ("VAR2") = s, Rcpp::Named("stringsAsFactors") = false ); } * NumericVector [#r8c3eeb8] - 旧来の場合 // 旧来の場合 SEXP mySum( SEXP x){ // 戻り値を用意 SEXP ans; int i, nx; // SEXPは直接操作できないので,C固有のポインタを用意 double *rx = REAL(x), *rans; nx = length(x); // ans を初期化.R側で解放されないようプロテクト PROTECT(ans = allocVector(REALSXP, 1)); // ans の操作は,そのポインタを保持するCのポインタを使う rans = REAL(ans); for(i = 0; i < nx; i++) { // Cのポインタを操作することで,ans を操作する rans[0] += rx[i]; } UNPROTECT(1);//プロテクトを外す return(ans); } # R CMD SHLIB mySum.cpp > dyn.load ("mySum.so") > x <- .Call("mySum", c(1.1, 1.2, 1.3, )) -- inline を使う場合 mySum <- cxxfunction (signature(x = "numeric") , ' int i, nx; double *rx = REAL(x), *rans; SEXP ans; nx = length(x); PROTECT(ans = allocVector(REALSXP, 1)); rans = REAL(ans); for(i = 0; i < nx; i++) { rans[0] += rx[i]; } UNPROTECT(1); return(ans); ' ) x <- mySum (c(1.1, 1.2, 1.3)) x - Rcppオリジナルのインターフェイスを使う場合 SEXP sum( SEXP x_ ){ Rcpp::NumericVector x(x_) ; double res = 0.0 ; for( int i=0; i<x.size(), i++){ res += x[i] ; } return Rcpp::wrap( res ) ; } - Rcpp sugar を使う SEXP sum( SEXP x_ ){ NumericVector x(x_) ; double res = sum( x ) ; return wrap( res ) ; } * CharacterVector [#pa7dd8c1] [[ここ:http://stackoverflow.com/questions/4106174/where-can-i-learn-to-how-to-write-c-code-to-speed-up-slow-r-functions]] にあるマンマだが,確かに手素が省ける.ただし mkChar よりも mkCharCE を使うべきかも知れない. - 旧来のインターフェイスを使う場合 SEXP foobar(){ SEXP ab; PROTECT(ab = allocVector(STRSXP, 2)); SET_STRING_ELT( ab, 0, mkChar("foo") ); SET_STRING_ELT( ab, 1, mkChar("bar") ); UNPROTECT(1); } - Rcppを使う.その1 RcppExport SEXP foobar(){ return Rcpp::CharacterVector::create( "foo", "bar" ) ; } - Rcppを使う.その2 RcppExport SEXP foobar(){ Rcpp::CharacterVector res(2) ; res[0] = "foo" ; res[1] = "bar" ; return res ; } - 以下でちょっと詰まった RcppExport SEXP testChar (SEXP a, SEXP b){ Rcpp::CharacterVector xa (a); Rcpp::CharacterVector xb (b); 続けて以下のようにしようとしたがコンパイラを通らない std::string sxa = xa[0]; 以下のようにしないと行けない std::string sxa; sxa = xa[0]; std::string sxb; sxb = xb[0]; * rcpp のサンプル [#xada24e9] Windows 7(64bit) 上の R-2.12.0 環境で,Rcpp_sample.cpp を以下のMakevars を用意. PKG_CXXFLAGS=$(shell Rscript -e "Rcpp:::CxxFlags()") PKG_LIBS=$(shell Rscript -e "Rcpp:::LdFlags()") してコマンドプロンプトから素直に R CMD SHLIB Rcpp_sample.cpp http://cygwin.com/cygwin-ug-net/using.html#using-pathnames g++ -shared -s -static-libgcc -o Rcpp_sample.dll tmp.def Rcpp_sample.cpp C:/Users/ishida/Documents/R/win-library/2.12/Rcpp/lib/i386/libRcpp.a -LC:/PROGRA~1/R/R-212~1.0/bin/i386 -lR とすると i386 が呼ばれる.確かに Rscript -e "Rcpp:::LdFlags() C:/Users/ishida/Documents/R/win-library/2.12/Rcpp/lib/i386/libRcpp.a となっている.が,こっちを呼んで欲しい Rscript --arch x64 -e "Rcpp:::LdFlags() C:/Users/ishida/Documents/R/win-library/2.12/Rcpp/lib/x64/libRcpp.a やむなく PKG_LIBS = C:/Users/ishida/Documents/R/win-library/2.12/Rcpp/lib/x64/libRcpp.a あるいは PKG_LIBS = -LC:/Users/ishida/Documents/R/win-library/2.12/Rcpp/lib/x64/ -lRcpp と修正し R --arch x64 CMD SHLIB Rcpp_sample.cpp と実行. dyn.load("Rcpp_sample.dll")