CでRを拡張するための方法解説

解説というほどの解説でもないんだけど、日本語リソースが非常に少ないので、ちょっとは役に立つかもしれないなと思って。

短かいですが、上のエントリのCの部分を中心に説明していこうかなと思います。

拡張するための3つ方法

Writing R Extensionsをばーっと読んで僕が把握したことを書きつらねていこうかと思います。4章以降ですね。

Rのシステムと他言語のインターフェイスとなる関数には大きくわけで3つ(4つ?)あります。

.C(.Fortran)
.Call
.External

の3つです。

Cの関数で書く

今回は自分が使った.Callのパターンについて書いてみようと思います。

R関数の実行を加速するためにCコードを使うことはしばしば非常に有効である。伝統的にはこれはRの.C関数を経由して行われる。このインタフェイスの一つの制限はRオブジェクトはC中で直接には処理できないことである。

とある通り、.CallはCの中でRのオブジェクトを扱いたい時に使います。僕の上のエントリのコードだと返り値をRの行列として返したかったので、.Cではなく.Callを使っています(そうじゃなくてもできそうだけどとりあえず)。RからこのCで書いた関数を呼び出す時には

dyn.load("/tmp/hoge.so")
r2norm <- function(n){
  data.frame(t(.Call("r2norm", as.integer(n))))
}

という風にやります。肝は

.Call("r2norm", as.integer(n))

の部分です。

具体的なCのコードの中身

.Callによる呼び出しをするときには関数の引数はSEXP型のオブジェクトとなっています。

処理すべき全てのRオブジェクトは型SEXPを持つ、これは定義型SEXPRECを持つ構造体へのポインタである。これはRオブジェクトの通常の全てのタイプ、つまり様々なタイプのベクトル、関数、環境、を処理できる可変な型と考えると分かりやすい。

とあって、難しいように思いますが、

詳細はこの節の後で与えられるが、多くの目的に取って、プログラマはそれを知っている必要はない。

とある通り、いくつか典型的な使い方を見てから理解するというのがいいかと思います。というか僕もあんまり分かっていないというw。

PROTECTとUNPROTECT

いくつか特徴的なところはあると思いますが、PROTECTの付近が一番特徴的かなと思うのでその辺を書いてみます。

もしCコード中にRオブジェクトを作ったら、オブジェクトへのポインタにPROTECTマクロを用いRにオブジェクトを使用することを知らせる必要がある。これはRにそのオブジェクトが使用中であることを告げ、そのためそれが破壊されることはなくなる。

と書いてあるようにCの中で(引数で受け取ったSEXP以外で)自分で作ったようなRオブジェクト(ここでは最終的にRで言うところのMatrixになっている)があれば、PROTECTしておきます。で、最後にPROTECTしているオブジェクトの数を引数にしてUNPROTECTしてあげないといけません。こっちが数えないといけないのはちょっと面倒ですが。。。

乱数発生時の注意点

Rは様々な乱数が生成できるので、ついCでも使いたくなるのですが、そういうことができます。今回の例ではRの正規乱数を生成してくれるものを使っています。下の部分です。

    prevX1 = 1.0 + 0.7 * (prevX2 - 2.0) + (norm_rand() * (1.0 - 0.7 * 0.7));

Rではrnorm(1,5,10)とやれば、平均が5標準偏差が10に従うような正規乱数が(この場合では)1つ生成できますが、上のnorm_randではそれらが指定できません。なので

mean + norm_rand() * var

という風にしてやればとりあえず目的のものは得られるんじゃないかなと思います。

乱数生成で注意しないといけない点は乱数生成前と生成後に書いてあげないといけないことがあって、それは下のようなものです。

  GetRNGstate();
  mean + norm_rand() * var;
  PutRNGstate();

これを挟んでやらないと毎回作らる乱数が同じになってしまうので注意しましょう。

ちなみにC側から呼び出せるR側の乱数生成の関数として

double unif_rand(); 
double norm_rand(); 
double exp_rand();

など(一様乱数、正規乱数、指数分布からの乱数)が与えられているようです。

キャスト関係

たぶんこの辺で一番時間がかかりました。引数にSEXP型の変数numがあるんですが、Rから呼び出されるときは10とかのただの数値です。がRのオブジェクトに一応なっているらしく明示的なキャストが必要となります。今回の例だとこの部分。

  int n = INTEGER(num)[0];

10とかだけなのでスカラーっぽい感じなんですが、Rには他言語の意味でのスカラーはありません。実際はベクトルになっています。なので[0]などとして要素を引っぱってきています。ちょっと嫌かな。。。

あと、Cで言うところのdoubleをRの(今回の例では)行列のオブジェクトansに代入していくとき、R側から見たdoubelにしてあげないといけません。そういうときのキャストはREALというのでやるようです。こんな感じ。

    REAL(ans)[i] = prevX1;

この辺はいくつかサンプルを見ながら、自分で使える例を貯めていくしかないですね。現在はWriting R Extensionsの例を見たりだとか、Rのソースを落としてきて、grepで使えそうなところを見つけたりだとかそんな感じです。