RのC拡張を触ってみたら修論用のプログラムが爆速になって新年早々鼻血がでそうになった件

以前からRにはC拡張というのがあって、そいつを利用すると凄まじい事になるという話を聞いていたのですが、色々あって後回しにしていました。
しかし、お正月なので（？）、ふと思い立って触ってみました。
意外とまとまった解説はなさそう（？）なので、今日は一番シンプルなインタフェイス関数.Cを利用する方法について紹介してみます。

参考

いきなりですが、参考にしたところ
http://d.hatena.ne.jp/syou6162/20090117/1232120983
Rから他言語利用 - RjpWiki
http://cran.r-project.org/doc/contrib/manuals-jp/R-exts.jp.pdf
というかWriting R Extensionsの4章を読めば何とかなる！！

手順

まずは大まかな手順を説明します。

RからCに投げたい部分の関数を作成する。
Cの関数をRに読み込む共用ライブラリとしてコンパイルする。
Rから共用ライブラリをロードする。
共用ライブラリに値を渡して処理してもらう。

以上です。
Rの知識というよりは、ちょいちょいポインタが絡んでくるので、C言語の多少の知識があった方が楽になる気がしました。

今回はRで行列を生成=>Cに渡して行列の要素全てに3を足す
というすごく適当な例を紹介してみたいと思います。

RからCに投げたい部分の関数を作成

いきなりコード

void mymatrix(double *a, int *ra, int *ca)
{
  int i, j;
  
  for (i = 0; i < *ra; i++) {
    for (j = 0; j < *ca; j++) {
      a[i*(*ca) + j] += 3;
    }
  }
}

これをmymatrix.cという名前で保存します。
ポイントは

関数の返り値をvoid型にすること
引数はポインタで与えること（RとCの変数の対応表がWriting R Extensionにあります。）
（今回の例では出てきませんが、）返り値にしたいものも引数に含めること

です。

今回の例では行列aと行列の行raと列caを引数にしています。

αとして行列は２次元配列のポインタと同じように先頭のポインタを渡して処理していく仕様になっているようです。

これはC言語で２次元配列を処理する時と同じ形になっているので、慣れている人には直感的に理解できるところですが、慣れてないと意外と詰まりそうなところですね。

Cの関数をRに読み込む共用ライブラリとしてコンパイル

環境はMacですが、

%R CMD SHLIB mymatrix.c

とすると何やらうにょうにょと処理してくれて、mymatrix.oとmymatrix.soというファイルができます。

Rから共用ライブラリをロードする。

これは簡単

dyn.load("./mymatrix.so")

と.soファイルをロードします。

共用ライブラリに値を渡して処理してもらう。

RのデータをCの関数に渡します。
ここでもいきなりコードを

dyn.load("mymatrix.so")

mymatrix <- function(a) {
  .C("mymatrix",
     result = as.double(a),
     as.integer(nrow(a)),
     as.integer(ncol(a)))$result
}
     
a <- matrix(1:10, 2, 5)

result <- matrix(mymatrix(a), 2, 5)

とすると、

> a
[,1] [,2] [,3] [,4] [,5]
[1,] 1 3 5 7 9
[2,] 2 4 6 8 10
> result
[,1] [,2] [,3] [,4] [,5]
[1,] 4 6 8 10 12
[2,] 5 7 9 11 13

となり、期待通りに動いていることが分かります。
ここでのポイントは、

.Cの第１引数はCの関数名とすること
第２引数以降は実際に関数に渡す値を与えること（但し、型のチェックを厳格に行う必要があります。）
$resultなどを利用しないと関数に与える引数全てをリストで返すことになるので、特定の値のみが欲しい場合は上の例のようにすればよいです。
ポインタとして渡しているので、結果はベクトル？一次元の行列？で返ってきます。必要に応じて元の行列の形に戻す必要があります。

結局何がおいしいか

ここに紹介した例では全くうまみはないです。
ただ処理が複雑になっただけw
ではどのような時にこのように一部の処理をCの関数に渡せばよいか？
それは本質的にfor文が要求されるような場面です。
２重３重のfor分が要求される場面では、Rは非常に遅いのでCに渡す事でパフォーマンスが劇的に改善します。

例えば、研究の中でグラフ上の全点の最短経路を求めたい場合に、Foyd-Warshall(Warshall-Foyd)法を利用しますが、これはfor文の3重ループが必要になります。
今までは速度が犠牲になることを覚悟の上で、他の処理の部分の楽さの方にメリットがあると判断して、Rを使用していましたが、これが実にボトルネックになっていました。
今回

void FW(double *data, int *NUM)
{
  int i, j, k;

  for (k = 0; k < *NUM; k++) {
    for (i = 0; i < *NUM; i++) {
      for (j = 0; j < *NUM; j++) {
        if (data[i*(*NUM) + j] > data[i*(*NUM) + k] + data[k*(*NUM) + j]) {
          data[i*(*NUM) + j] = data[i*(*NUM) + k] + data[k*(*NUM) + j];
        }
      }
    }
  }
}

というC拡張を書き、

FoydWarshall <- function(data) {
  dyn.load("./FW.so")
  CFoydWarshall <- function(data) {
    print("Culculate The Shortest Path")
    .C("FW",
       result = as.double(data),
       as.integer(NUM))$result
  }
  result <- CFoydWarshall(data)
  result[which(result == 1000)] <- 0
  return(result)
}

とRからCに投げるようにしたところ、パフォーマンスがあり得ないくらい向上しました。
これはすごい。マジで昨日の深夜に鼻血がでそうなくらい興奮してしまいました。

まとめ

最近意識的にRで研究用のプログラムを書いていますが、時々妙なところで躓くので、自分的なまとめとして、

まずRの組み込み関数を探す（パッケージがあるかもしれないのでやはり探す）
上記では微妙に手が届かない場合は、行列演算や、ベクトル演算に持ち込むように処理を工夫する
どうしてもfor文でしか実現できなさそうな処理は、まずそこがボトルネックになっているかどうかを判断。ボトルネックになっていそうならCに投げる

というのが速度面に関しては大事な心構えだと思います。

おまけ

しかし、上記のインタフェイス関数.Cでは実に単純な処理しかできない。
例えば上の例で

result[which(result == 1000)] <- 0

としている部分で、本来はInfとしたいところですが、Rの無限大を表すInfをCでは使えないため、やむなく大きな数として1000を代入しました。
こういう風にRでは普通にできることがCになるとできなくなって、若干不便だなーということが出てきます。
そこで、CでRのデータ構造（行列とかベクトルとか欠損値とかなんやらかんやら）を扱うために、.Call関数と.External関数が用意されていたりするわけです。
その場合、先ほどは関数しかなかった.Cファイルに

#include<R.h>
#include<Rinternals.h>

というヘッダーをインクルードします。
Rinternals.hを見れば分かるのですが、様々なマクロが定義されており、それを通じてRのデータ構造を扱う事ができるようになるようです。
先ほどの「正の無限大」は"R_PosInf"のようです。

というように、現在勉強中だけど、こいつもなかなか面白い感じです。
またまとまった知識がついたらここで紹介するかもです。