トップ 差分 一覧 ソース 検索 ヘルプ RSS ログイン

BugTrack-R備忘録/56

R備忘録 /状態空間モデリング/donlp2/その他のメモ

R備忘録 - 記事一覧

R(R言語)からCを呼ぶ(Linux編)

  • 投稿者: みゅ
  • カテゴリ: なし
  • 優先度: 普通
  • 状態: 完了
  • 日時: 2010年01月23日 19時40分27秒

内容

  • RからCを呼ぶ
    • Rはなにかアイデアとか思いついたことをすぐに試すにはよいし、とても使いやすいのですが、さすがにインタープリターだけあって、C言語なんかに比べるととても遅いのです.そんなときは時間のかかる部分をCとかで書いてRからそれを呼び出すということができます.はい、ばりばり使いまくってます.
  • とにかくめちゃくちゃ速いです.いやRが遅いというべきか?
  • バグ報告や質問は別館でお願いします.

Makefileを用意する

  • Makefileを用意します.以下のリンクなんかも参考にしてみてください.
    • BugTrack-R備忘録/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-R備忘録/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
//R用 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
  • 1番目の引数がC関数の名前で、その後はCで作成した関数の引数の順番どおりに引数を記述します.
> conv(1:3,2:6)
[1]  2  7 16 22 28 27 18
  • となり、きちんとCが呼べています.

サンプル

  • BugTrack-R備忘録/53の5.9のサンプルです.今度は「.Call」関数です.
  • 先ほどのtest.cにつづけて書いてしまってかまいません.
  • <R.h>と<Rdefines.h>も追加してください.
//R用 C
#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」を使うことが多いですね.

コメント