Rcpp - アールメカブ

アールメカブ


Rcpp

Rの備忘録

dirk.eddelbuettel Rcpp

_ 基本

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

Rcpp::NumericVector xa(a);

この xa はそのまま演算できる.

  1. RcppExport SEXP convolve3cpp(SEXP a, SEXP b) {
  2.   Rcpp::NumericVector xa(a);
  3.   Rcpp::NumericVector xb(b);
  4.   int n_xa = xa.size(), n_xb = xb.size();
  5.   int nab = n_xa + n_xb - 1;
  6.   Rcpp::NumericVector xab(nab);
  7.   for (int i = 0; i < n_xa; i++)
  8.     for (int j = 0; j < n_xb; j++)
  9.       xab[i + j] += xa[i] * xb[j];
  10.   return xab;
  11. }

これをRに返す際には,明示的に型変換する必要がないので,ついつい SEXP そのものだと誤解しそうになるが,あくまで Rcpp::RObjectクラスのオブジェクトであることを忘れないようにしないと,妙な代入などをおこなってしまいそうになる.

ちなみに,伝統的な R API で,Rinternals.h を使う場合は,演算処理のためにポインタをとる.

  1. SEXP convolve2(SEXP a, SEXP b)
  2.      {
  3.          R_len_t i, j, na, nb, nab;
  4.          double *xa, *xb, *xab;
  5.          SEXP ab;
  6.           // 型を指定して保護
  7.          PROTECT(a = coerceVector(a, REALSXP));
  8.          PROTECT(b = coerceVector(b, REALSXP));
  9.          na = length(a); nb = length(b); nab = na + nb - 1;
  10.          PROTECT(ab = allocVector(REALSXP, nab));
  11.     // C言語での演算のためポインタを取得
  12.          xa = REAL(a); xb = REAL(b);
  13.          xab = REAL(ab);
  14.          for(i = 0; i < nab; i++) xab[i] = 0.0;
  15.          for(i = 0; i < na; i++)
  16.              for(j = 0; j < nb; j++) xab[i + j] += xa[i] * xb[j];
  17.          UNPROTECT(3);
  18.     //SEXP型なのでそのまま返却できる
  19.          return(ab);
  20.      }

Rdefines.h を使う場合も,演算処理のためにポインタをとる.

  1. SEXP convolve2(SEXP a, SEXP b)
  2.      {
  3.          R_len_t i, j, na, nb, nab;
  4.          double *xa, *xb, *xab;
  5.          SEXP ab;
  6.  
  7.     // 数値型として保護し
  8.          PROTECT(a = AS_NUMERIC(a));
  9.          PROTECT(b = AS_NUMERIC(b));
  10.          na = LENGTH(a); nb = LENGTH(b); nab = na + nb - 1;
  11.          PROTECT(ab = NEW_NUMERIC(nab));
  12.          // Cの演算のためにポインタを取り出す
  13.          xa = NUMERIC_POINTER(a); xb = NUMERIC_POINTER(b);
  14.          xab = NUMERIC_POINTER(ab);
  15.          // ポインタを通した計算
  16.          for(i = 0; i < nab; i++) xab[i] = 0.0;
  17.          for(i = 0; i < na; i++)
  18.               for(j = 0; j < nb; j++) xab[i + j] += xa[i] * xb[j];
  19.          UNPROTECT(3);
  20.          // SEXP型なのでそのまま返却できる
  21.          return(ab);
  22.      }

と明示的にPROTECTが必要になる.

_ 引数を受ける

こんな感じで利用するのが一番楽?

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? 作成

#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?

  • 旧来の場合
    // 旧来の場合
    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?

ここ にあるマンマだが,確かに手素が省ける.ただし 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 のサンプル

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")
 
Link: R_old_tips3(1744d) Rの備忘録(3781d)
Last-modified: 2011-09-14 (水) 17:05:44 (4607d)