2015-01-13

Rcppおさらい

R Rcpp Cpp

Rcppを使ってみた(こちら)
おおまかに分かったので、じゃあ、Rcpp/Cppの強みについておさらいしておこう
- 一番の基本(この記事)
- 少し先(sugar,STL,Armadillo)
こちらを読む
Rcppをいつ使う？
- Rのコードでできるだけのことをしても「遅い」とき
- それはいつ？
  - ループをRが得意なベクトル化で高速処理できないとき(前の結果に依存するループとか)
  - 再帰関数や、関数呼び出しが繰り返しなされるとき(関数呼び出しのオーバーヘッドはC++ << R)
  - データ構造としてC++では使いやすくて、Rでは使いにくいものを扱うとき(standard template library (STL)というのは、C++が持つデータ構造ライブラリ→参考１,参考２,参考３)
Rcppパッケージは、(主に)２つの方法でC++関数をコンパイルしてRから呼び出して使えるようにする
- 短いならcppFunction()で可、長ければsourceCpp()を使うのが普通だろう。どちらか片方を覚えるならsourceCpp()だろう。
- (1) cppFunction()関数を使う
  - C++関数を文字列としてcppFunction()に渡す
- (2) sourceCpp()関数を使う
  - C++関数をファイルで書いてsourceCpp()関数に読み込ませる
- ２整数の積を返す関数
  - R

rprod <- function(x,y){
  x*y
}
x <- 3
y <- 2
rprod(x,y)

- - - 実行結果

> rprod <- function(x,y){
+     x*y
+ }
> x <- 3
> y <- 2
> rprod(x,y)
[1] 6

- - cppFunction()を使う

prodstring <- 'int add(int x, int y) {
  int pr = x * y;
  return pr;
}'
cppFprod <- cppFunction(prodstring)
cppFprod
cppFprod(x,y)

- - - 実行結果

> prodstring <- 'int add(int x, int y) {
+   int pr = x * y;
+   return pr;
+ }'
> cppFprod <- cppFunction(prodstring)
> cppFprod
function (x, y) 
.Primitive(".Call")(<pointer: 0x0000000064541770>, x, y)
> cppFprod(x,y)
[1] 6

- - sourceCpp()を使う
    - cpp関数のファイルをprcpp.cppとして作成保存する
      - cppFunction()関数は、このc++ファイルのヘッダー部分や//Rcpp::exportなどをやってくれるように書かれている関数とわかる

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
int prcpp(int x, int y) {
  int pr = x * y;
  return pr;
}

- - - 読み込み実行

sourceCpp("prcpp.cpp")
prcpp
prcpp(x,y)

- - - 実行結果

> sourceCpp("prcpp.cpp")
> prcpp
function (x, y) 
.Primitive(".Call")(<pointer: 0x00000000629c1770>, x, y)
> prcpp(x,y)
[1] 6

R関数とc++関数の違い
- 両者を並べてみる

rprod <- function(x,y){
  x*y
}

int prcpp(int x, int y) {
  int pr = x * y;
  return pr;
}

- Rでは関数を作成する関数function()を使うが、C++では使わない
- Rではreturn()を使っても使わなくてもよいが、C++では必ず使う
- 行末記号をRは使わないがC++は使う
- 型宣言するところが違う
  - 関数prcpp()の帰り値がintであることを"int prcpp(..."のように宣言
  - 引数がintであると宣言(int x, int yのようにx,yがintであると)
  - 関数内部でも新出の変数の型を宣言 int pr ...
  - これはC++仕様の基本なのでどこまでもついて回る
速さ比べ
- 格子酔歩１
  - 原点から出発して二次元格子の４方向を等確率で選んで動く酔歩、n歩分のシミュレーション
  - Rの場合

n <- 10000
rRW <- function(n){
  stepx <- sample(c(-1,1),n,replace=TRUE)
  stepy <- sample(c(-1,1),n,replace=TRUE)
  xORy <- sample(c(0,1),n,replace=TRUE)
  xselect <- which(xORy==0)
  yselect <- which(xORy==1)
  x <- y <- rep(0,n)
  x[xselect] <- stepx[xselect]
  y[yselect] <- stepy[yselect]
  x <- c(0,x)
  y <- c(0,y) 
  return(c(cumsum(x),cumsum(y)))
}
outR <- rRW(n)
plot(matrix(outR,ncol=2),type="l")

- - c++の場合

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector cppRW(int n){
  NumericVector ret((n+1)*2);
  for(int i=1;i<(n+1);i++){
    int tmp = n + i + 1;
    ret[i] = ret[i-1];
    ret[tmp] = ret[tmp-1];
    if(rnorm(1)[0] < 0){
      if(rnorm(1)[0] < 0){
        ret[i]=ret[i]+1;
      }else{
        ret[i]=ret[i]-1;
      }
    }else{
      
      if(rnorm(1)[0] < 0){
        ret[tmp]=ret[tmp]+1;
      }else{
        ret[tmp]=ret[tmp]-1;
      }
    }
  }
  return ret;
}

- - - 実行

sourceCpp("cppRW.cpp")
n <- 10000
outCpp <- cppRW(n)
plot(matrix(outCpp,ncol=2),type="l")

- 比較

system.time(cppRW(n*100))
system.time(rRW(n*100))

- - 結果。大差なし

> system.time(cppRW(n*100))
   ユーザ   システム       経過  
      0.31       0.02       0.33 
> system.time(rRW(n*100))
   ユーザ   システム       経過  
      0.34       0.02       0.36

- 酔歩は「過去」に依存して発展するとは言え、４方向のどちらに進むかについて、過去に影響されずに選べるので、ステップ方向の乱択がベクトル化で来ている
- では、Self-avoiding pathのようにステップ方向を決めるときに、過去の経路を確認しないといけなくするとどうなるか、やってみる。ループを途中で切り上げたりとか工夫もできるが、べたにやる。
  - Rバージョン

rRWsa <- function(n){
  stepx <- sample(c(-1,1),n,replace=TRUE)
  stepy <- sample(c(-1,1),n,replace=TRUE)
  xORy <- sample(c(0,1),n,replace=TRUE)
  xselect <- which(xORy==0)
  yselect <- which(xORy==1)
  x <- y <- rep(0,n)
  x[xselect] <- stepx[xselect]
  y[yselect] <- stepy[yselect]
  X <- rep(0,n+1)
  Y <- rep(0,n+1)
  for(i in 2:(n+1)){
    tmpx <- X[i-1]
    tmpy <- Y[i-1]
    tmpx2 <- tmpx + x[i-1]
    tmpy2 <- tmpy + y[i-1]
    tmpx3 <- X[1:(i-1)]-tmpx2
    tmpy3 <- Y[1:(i-1)]-tmpy2
    tmpD <- abs(tmpx3) + abs(tmpy3)
    if(min(tmpD)>0){ # 条件に合わなければ進まない
      X[i] <- tmpx2
      Y[i] <- tmpy2
    }else{
      X[i] <- X[i-1]
      Y[i] <- Y[i-1]
}
  }
  return(c(X,Y))
}
n <- 10000
outRsa <- rRWsa(n)
plot(matrix(outRsa,ncol=2),type="l")

- - CPP

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector cppRWsa(int n){
  NumericVector ret((n+1)*2);
  for(int i=1;i<(n+1);i++){
    int tmp = n + i + 1;
    ret[i] = ret[i-1];
    ret[tmp] = ret[tmp-1];
    int tmpx = ret[i];
    int tmpy = ret[tmp];
    if(rnorm(1)[0] < 0){
      if(rnorm(1)[0] < 0){
        tmpx=tmpx+1;
      }else{
        tmpx=tmpx-1;
      }
    }else{
      
      if(rnorm(1)[0] < 0){
        tmpy=tmpy+1;
      }else{
        tmpy=tmpy-1;
      }
    }
    bool ok = TRUE;
    for(int j=0;j<i;j++){
      if(abs(ret[j]-tmpx)+abs(ret[n+1+j]-tmpy) == 0){
        ok = FALSE;
      }
    }
    if(ok){
      ret[i] = tmpx;
      ret[tmp] = tmpy;
    }
  }
  return ret;
}

- - 比較

sourceCpp("prcpp.cpp")
n <- 100000
system.time(cppRWsa(n))
system.time(rRWsa(n))

- - 結果

> n <- 10000
> 
> system.time(cppRWsa(n))
   ユーザ   システム       経過  
      0.09       0.00       0.09 
> system.time(rRWsa(n))
   ユーザ   システム       経過  
      0.92       0.00       0.92 
> 
> n <- 100000
> 
> system.time(cppRWsa(n))
   ユーザ   システム       経過  
     10.31       0.00      10.33 
> system.time(rRWsa(n))
   ユーザ   システム       経過  
     96.39       0.29      97.59

Listの出し入れ
- ベクトル・行列の出し入れは比較的簡単なはずで、確かにその通り
- Rのリストはいろんなデータ型が入り乱れているので、それに対応してくれるようにRcppは設計してある

// [[Rcpp::export]]
double hoge(List mod) {
...

Rの関数(f)をC++に渡すこともできる
- ここでcpp関数の返り値のデータ型が"RObject"となっている。Rの関数f()がどんなものかわからないので、何にでも対応(Listが一番自由度が高い)するようになっている、ということ(らしい)

// [[Rcpp::export]]
RObject hage(Function f) {

基本は終わり
- 次の記事は、便利に使う、楽に使う、より上手に使う、ためにsugar,STL,Armadilloなどについて