2012-01-03

RcppをWindowsで使う

RcppはRの(関数を作ったりするときに使う)パッケージでC++を使うためのパッケージ
WindowsでRを使っており、ちょっと前作業も必要なので、冬休みくらいに、時間があるときに向いている
こちらで書かれているように「むっちゃ速い」ということなので、試してみよう
むっちゃ速いというサイトのソースは、コピーペースト-readyではないので、その本家のサイト(こちら)のソースを使ってみる
「Cってどんな感じ？」と試したのが数年前(３日間とか？）、C++っていいよ、と言われて試したこともないけど、その程度の知識でも使えるようになるか、やってみよう
その他参考(こちらやこちら→と思ったら、後者は読めないサイトでした…)

# Xian's code, using <- for assignments and passing x down
# 単純な計算をn回実施する関数５種
f <- function(n, x=1) for (i in 1:n) x=1/(1+x)
g <- function(n, x=1) for (i in 1:n) x=(1/(1+x))
h <- function(n, x=1) for (i in 1:n) x=(1+x)^(-1)
j <- function(n, x=1) for (i in 1:n) x={1/{1+x}}
k <- function(n, x=1) for (i in 1:n) x=1/{1+x}

# now load some tools
# インストールして読み込む
# ２種類使うらしい
# 基本のパッケージは"Rcpp"。Rcppのソースをコンパイルして…という処理をしてくれる諸処理のパッケージ
# inlineパッケージはRのコマンドでC,C++ソース(のようなもの)を書いてC,C++関数オブジェクトにする関数を持つ
# "Rcpp"をプラグインとして使って、そうするらしい
install.packages("Rcpp")
install.packages("inline")
library(Rcpp)
library(inline)

# and define our version in C++
# inline パッケージのcxxfunction()関数を使ってC++関数オブジェクトを作る
# signature(ns="integer",xs="numeric")は関数の変数宣言
# 'int ... wrap(x); 'のようにC++のソース(多分)
# int n = as<int>(ns); double x=as<double>(x);で変数宣言して
# for (int ...以下に処理を書き
# return wrap(x)で返す
# plugin="Rcpp"でRcppを使ってコンパイル(等)をすることを示している(ようだ)
l <- cxxfunction(signature(ns="integer", xs="numeric"),
                 'int n = as<int>(ns); double x=as<double>(xs);
                  for (int i=0; i<n; i++) x=1/(1+x);
                  return wrap(x); ',
                 plugin="Rcpp")

これをコピーペーストしただけでは、次のようなエラーが出る

g++: not found
g++: not found CRAN windows

"コンパイルできない"まずは、コンパイル問題を解決しよう。gccがないとのこと、次のようなキーワードでウェブ検索するとこんなサイトに行きつく
- Rの64bit版と32bit版とでコンパイラに違いがあって、それに伴うトラブルもあるようだ(64ビット版のR x64-2.14.1ではだめだったがR 2.14.1 (32bit版)ではうまく行った)
- 以下のように解決する方法もあるが、こちらに書いたように、Rtoolsを使って解決することもできる
その中で必要とされているMinGWとMSYSで検索すると、それぞれが、コンパイラなどの一括ダウンロードセットのようなものであることがわかる
こちらを参考にしてインストールしてみる
- WinGWもMSYSも自動インストーラなので、問題ない
- コマンドプロンプトを立ち上げて、コマンド "gcc" を打つと、「動かない」
- WindowsのPATHを変えることが必要なので、こちらを参考にしてPATHを変える

- こんな画面が出て、ハイライトしているところを書き換えるわけだけれど、「初心に戻って」考えると、こういうときに悩むのは、『本当に書き換えても大丈夫なのか』ということ
- というわけで、エクスプローラで拡張子を表示させるという設定変更にも抵抗がある人向けの説明を書こう
- まず、ハイライト部分をすべてハイライトしてあることを確認した上で、コピーしてテキストエディタにペーストする
- 中身がよく観察できる環境にした上で、「さて、パスを通すとは、どうすればいいか」をよく眺める。良く眺める理由は、いざとなったら、自分の責任でもとに戻せるため
- ';'で継ぎ足し、継ぎ足しをした構成になっていることを自分の眼で確認したら、あとは、継ぎ足すべき文字列を決める
  - さきほどインストールしたMinGWは"C:\MinGW"に入ったのですが、それを確認した上で、そこの"bin"フォルダを確認しつつ、その中も見ておきます

- - このフォルダをしていすればよいので、次のように書き換わります

...;C:\MinGW\bin;C:\msys\1.0\bin

こうしたうえで、コマンドプロンプトを立ち上げなおして、"gcc"とかのコマンドを打つと、今度は、(コンパイルされるべきファイルを渡していないので、エラーメッセージは出ますが)"gcc"というコマンド自体は認識されていることがわかります
さて、念のためRを再起動して、もう一度、以下を回すと、今度は動きます

l <- cxxfunction(signature(ns="integer", xs="numeric"),
                 'int n = as<int>(ns); double x=as<double>(xs);
                  for (int i=0; i<n; i++) x=1/(1+x);
                  return wrap(x); ',
                 plugin="Rcpp")

この処理でできるオブジェクトの中身を見てみよう

> l
An object of class "CFunc"
function (ns, xs) 
.Primitive(".Call")(<pointer: 0x66781290>, ns, xs)
<environment: 0x09809ec8>
Slot "code":
[1] "\n// includes from the plugin\n\n#include <Rcpp.h>\n\n\n#ifndef BEGIN_RCPP\n#define BEGIN_RCPP\n#endif\n\n#ifndef END_RCPP\n#define END_RCPP\n#endif\n\nusing namespace Rcpp;\n\n\n// user includes\n\n\n// declarations\nextern \"C\" {\nSEXP filea6034c539e9( SEXP ns, SEXP xs) ;\n}\n\n// definition\n\nSEXP filea6034c539e9( SEXP ns, SEXP xs ){\nBEGIN_RCPP\nint n = as<int>(ns); double x=as<double>(xs);\n                  for (int i=0; i<n; i++) x=1/(1+x);\n                  return wrap(x); \nEND_RCPP\n}\n\n\n"

Rcppパッケージ、inlineパッケージはうまく動いたようなので、次に進みます

# more tools
# 何かしらパッケージを使っている
library(rbenchmark)

# now run the benchmark
# Rcppの「威力」を試すための便利パッケージのように見える
# 「大きい数」を与えて、「たくさんの処理」をさせる
N <- 1e6
# その「たくさんの処理」の実際
# Cpp化していない関数f(),g(),h(),j(),k()とCpp化した関数l()とを並べている処理なので
# Cpp化した関数が「速い」ことを示す処理のはず
benchmark(f(N, 1), g(N, 1), h(N, 1), j(N, 1), k(N, 1), l(N, 1),
          columns=c("test", "replications", "elapsed", "relative"),
          order="relative", replications=10)

benchmark()の出力を見ると、何をする関数なのかはわかる
- 複数の関数に繰り返し処理をさせて、所要時間の測定と、相対的な速さを「速い順」に出す関数だ

> benchmark(f(N, 1), g(N, 1), h(N, 1), j(N, 1), k(N, 1), l(N, 1),
+           columns=c("test", "replications", "elapsed", "relative"),
+           order="relative", replications=10)
     test replications elapsed relative
6 l(N, 1)           10    0.29  1.00000
1 f(N, 1)           10    9.22 31.79310
4 j(N, 1)           10   10.73 37.00000
2 g(N, 1)           10   11.11 38.31034
3 h(N, 1)           10   13.34 46.00000
5 k(N, 1)           10   22.73 78.37931

思ったより手早く、終わった。年末の雑事をしながら、実質２時間の作業か。これなら、タイムパフォーマンスは悪くない
後は、次の部分を書くことが億劫でなくなるだけ

l <- cxxfunction(signature(ns="integer", xs="numeric"),
                 'int n = as<int>(ns); double x=as<double>(xs);
                  for (int i=0; i<n; i++) x=1/(1+x);
                  return wrap(x); ',
                 plugin="Rcpp")