//このページが表示される方は、URLから「action=SOURCE&」を削除してみてください
[[Ｒ備忘録 - 記事一覧]]

!!!R（R言語）からCを呼ぶ（Linux編）
*投稿者： みゅ
*カテゴリ： なし
*優先度： 普通
*状態： 完了
*日時： 2010年01月23日 19時40分27秒
//{{bugstate}}
!!内容
*RからCを呼ぶ
**Ｒはなにかアイデアとか思いついたことをすぐに試すにはよいし、とても使いやすいのですが、さすがにインタープリターだけあって、Ｃ言語なんかに比べるととても遅いのです．そんなときは時間のかかる部分をＣとかで書いてＲからそれを呼び出すということができます．はい、ばりばり使いまくってます．
*とにかくめちゃくちゃ速いです．いやRが遅いというべきか？
*バグ報告や質問は[別館|http://phoenixx.sakura.ne.jp/mt/R/2010/01/rc.html]でお願いします．
!!Makefileを用意する
*Makefileを用意します．以下のリンクなんかも参考にしてみてください．
**[[BugTrack-Ｒ備忘録/28]] - Rでロードできるライブラリ（.soファイル）のコンパイルのお作法？（Linux）
*Makefileサンプル
 R_SHARE_DIR = ${R_HOME}/share/
 R_INCLUDE_DIR = ${R_HOME}/include
 include ${R_HOME}/etc${R_ARCH}/Makeconf
 HEADER = test.h
 
 %.o: %.c $(HEADER)
 	$(CC) $(ALL_CPPFLAGS) $(ALL_CFLAGS) -c $< -o $@
 
 %.so : %.o
 	$(SHLIB_LINK) -o $@ $< $(ALL_LIBS) $(BLAS_LIBS) $(LAPACK_LIBS)
 
 clean :
 	rm -f *.o

!!サンプル①
*「[[BugTrack-Ｒ備忘録/53]] - Writing R Extensions(Version 2.7.1 (2008-06-23)):::5 システムと外部言語へのインターフェース」
*このあたりのサンプルから始めてみましょう．
 void convolve(double *a, int *na, double *b, int *nb, double *ab)
 {
     int i, j, nab = *na + *nb - 1;
     for(i = 0; i < nab; i++)
         ab[i] = 0.0;
     for(i = 0; i < *na; i++)
         for(j = 0; j < *nb; j++)
             ab[i + j] += a[i] * b[j];
 }
*なにをやっている関数なんでしょうね（笑．convolveだから畳み込み演算なんでしょうけど（笑
*とにかくこれでtest.cというファイルを作成します．
*test.c
 //Ｒ用 Ｃ
 
 void convolve(double *a, int *na, double *b, int *nb, double *ab)
 {
     int i, j, nab = *na + *nb - 1;
     for(i = 0; i < nab; i++)
         ab[i] = 0.0;
     for(i = 0; i < *na; i++)
         for(j = 0; j < *nb; j++)
             ab[i + j] += a[i] * b[j];
 }
*コンパイルして、その後Rから呼ぶことができるダイナミックリンクライブラリを作成するのですが、さっき作ったMakefileで、コマンド一発です．
 make test.so
*これによって、コンパイル>>>.soファイル作成全部やってくれます．らくちん．
 $ make test.so
 gcc -std=gnu99 -I/usr/local/lib/R/include  -I/usr/local/include    -fpic  -O2 -m
 tune=i686 -c test.c -o test.o
 gcc -std=gnu99 -shared -L/usr/local/lib -o test.so test.o   -L/usr/local/lib/R/l
 ib -lR -L/usr/lib/lapack/atlas -llapack -L/usr/lib/blas/atlas -lblas -L/usr/lib
 -latlas
 rm test.o
*さて出来上がった「.so」ファイルを呼びます．まずRを立ち上げましょう．
 dyn.load("test.so")
*「.so」があるディレクトリと同じ場所でRを立ち上げると、上のように打つだけで、「.so」が読み込まれています．ちゃんと関数の名前が認識できているか確認してみます．
 > is.loaded("convolve")
 [1] TRUE
*convolveというCのエントリを認識しています．
*さてこのCで作成した関数を呼ぶために「.C」というRの関数を使います．
 conv <- function(a, b)
      .C("convolve",
     as.double(a),
     as.integer(length(a)),
     as.double(b),
     as.integer(length(b)),
     ab = double(length(a) + length(b) - 1))$ab
*１番目の引数がC関数の名前で、その後はCで作成した関数の引数の順番どおりに引数を記述します．
 > conv(1:3,2:6)
 [1]  2  7 16 22 28 27 18
*となり、きちんとCが呼べています．
!!サンプル②
*[[BugTrack-Ｒ備忘録/53]]の５．９のサンプルです．今度は「.Call」関数です．
*先ほどのtest.cにつづけて書いてしまってかまいません．
*<R.h>と<Rdefines.h>も追加してください．
 //Ｒ用 Ｃ
 #include <R.h>
 #include <Rdefines.h>
 
 void convolve(double *a, int *na, double *b, int *nb, double *ab)
 {
     int i, j, nab = *na + *nb - 1;
     for(i = 0; i < nab; i++)
         ab[i] = 0.0;
     for(i = 0; i < *na; i++)
         for(j = 0; j < *nb; j++)
             ab[i + j] += a[i] * b[j];
 }
 
 SEXP convolve2(SEXP a, SEXP b)
 {
   int i, j, na, nb, nab;
   double *xa, *xb, *xab;
   SEXP ab;
   PROTECT(a = AS_NUMERIC(a));
   PROTECT(b = AS_NUMERIC(b));
   na = LENGTH(a); nb = LENGTH(b); nab = na + nb - 1;
   PROTECT(ab = NEW_NUMERIC(nab));
   xa = NUMERIC_POINTER(a); xb = NUMERIC_POINTER(b);
   xab = NUMERIC_POINTER(ab);
   for(i = 0; i < nab; i++) xab[i] = 0.0;
   for(i = 0; i < na; i++)
     for(j = 0; j < nb; j++) xab[i + j] += xa[i] * xb[j];
       UNPROTECT(3);
   return(ab);
 }
*このように書いてmake．Rを立ち上げて、
 > dyn.load("test.so")
 > is.loaded("convolve")
 [1] TRUE
 > is.loaded("convolve2")
 [1] TRUE
 > conv2 <- function(a, b) .Call("convolve2", a, b)
 > conv2(1:3,2:6)
 [1]  2  7 16 22 28 27 18
*はい、できあがり．
*こちらのタイプはRのオブジェクトをそのままCに渡すことができます．そのかわりCの中で、そのRのオブジェクトがきちんと扱えるようなマクロを駆使しなければなりませんが（苦笑
*どちらを使うかは、ケースバイケースでしょう．ですが管理人みゅは最近は「.Call」を使うことが多いですね．


!!コメント
//{{comment}}