Rcpp のバックアップの現在との差分(No.5) - アールメカブ

アールメカブ


Rcpp のバックアップの現在との差分(No.5)


  • 追加された行はこの色です。
  • 削除された行はこの色です。
[[Rの備忘録]]

[[dirk.eddelbuettel Rcpp:http://dirk.eddelbuettel.com/code/rcpp.html]]

#contents

* 基本 [#hbf6ee84]
Rcppインターフェイスを使ったクラスライブラリを作成するコード側で,たとえば
Rから受け取った SEXP 型 a は,以下のように Rcpp::RObjectクラスの派生クラスに変換する必要があるのだが,

 Rcpp::NumericVector xa(a);

この xa はそのまま演算できる.
#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 を使う場合は,演算処理のためにポインタをとる.

#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 を使う場合も,演算処理のためにポインタをとる.

#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){
   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++) {
 
  // 戻り値を用意
  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);
  }
  UNPROTECT(1);//プロテクトを外す
  return(ans);
 }

 x <- .Call("mySum.so", c(1.1, 1.2, 1.3, ))
 # 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]] にあるマンマだが,確かに手素が省ける
[[ここ: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 

 SEXP foobar(){
 RcppExport SEXP foobar(){
   return Rcpp::CharacterVector::create( "foo", "bar" ) ;
 }

- Rcppを使う.その2

 SEXP foobar(){
 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")