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

BugTrack-R備忘録/57

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

R備忘録 - 記事一覧

R(R言語)からCを呼ぶ(Linux編) - ちょっとしたコツ

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

内容

  • 「ちょっとしたコツ」集
  • 管理人みゅの備忘録(笑
  • ここから先は、BugTrack-R備忘録/56のCのソースにどんどん追加していくだけで、動くはずです.
  • ソース → test.c(13)
  • バグ報告や質問は別館でお願いします.

Cに渡したRオブジェクトのディメンジョンを得る

  • GET_DIM
  • 行列なんかを渡すときに、そのディメンジョンをいちいち整数で渡さなくてよい.Cの中で取得できる.
  • Cサンプル
SEXP test001(SEXP a){
  int n, m;
  n = INTEGER(GET_DIM(a))[0];
  m = INTEGER(GET_DIM(a))[1];
  printf("nrow of a : %d\n", n);
  printf("ncol of a : %d\n", m);
  return a;
}
  • Rから呼ぶ
> a <- matrix(double(1), 3,5)
> .Call("test001", a)
nrow of a : 3
ncol of a : 5
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    0    0    0    0    0
[3,]    0    0    0    0    0

dgesv(lapack)を使う

  • Solves a general system of linear equations AX=B.
  • http://www.netlib.org/lapack/double/dgesv.f
  • こんな計算はRでちょちょいのちょいなんですが、スピードアップを狙ってCで記述したいときに、こんな計算をするためにいちいちRに戻って計算させていては、それこそパフォーマンスが出ないので・・・.こんな計算でもCの中でやらせたほうがよいです.ですので、blasとかlapackは必須だったりします.
  • C言語
SEXP Rdgesv01(SEXP A, SEXP b, SEXP info)
{
  int i;
  long int n = INTEGER(GET_DIM(A))[0];
  long int nrhs=1;
  long int piv[n];   // こんな書き方できたっけ???

  printf("N = %ld\n", n);
  dgesv_(&n, &nrhs, REAL(A), &n, piv, REAL(b), &n, INTEGER(info));
  for(i=0; i<n; ++i) printf("%lf\n", REAL(b)[i]);
  return(b);
}
  • Rから、Ax=bを解く
A <- cbind(c(1, 3, 1), c(1, 1, -2), c(1, -3, -5))
b <- c(3,1,-6)
solve(A, b)
> solve(A, b)
[1] 1 1 1   # Rによる解
> info <- integer(1)
> .Call("Rdgesv01", A, b, info)
N = 3   #Cの中で表示させている
1.000000   #Cの中で表示させている
1.000000   #Cの中で表示させている
1.000000   #Cの中で表示させている
[1] 1 1 1   #「.Call」の戻り値
  • ただし今のCの書き方だと、引数にした「b」自体の中身も書き換わってしまっています.dgesvの仕様なので仕方ありません.(ついでに言うとAも書き換わる・・・)
> b
[1] 1 1 1
  • これを、防ぐには結果を入れる場所を用意して、そこに結果をセットして返すしかありません.
  • Rでこれを実現するのは簡単です.「A」「b」が既にセットされていると仮定すると
tmpA <- A
tmpb <- b
.Call("Rdgesv01", tmpA, tmpb, info)
  • 別の入れ物を作ってそれで計算させればいいだけです・・・
  • でもこれだと、Rで代わりの入れ物を用意しているので、Cの中で一連の処理をさせたい場合には使えません.そこでこれをCの中でちゃんと用意するようにします.
SEXP Rdgesv01(SEXP A, SEXP b, SEXP info)
{
  int i;
  long int n = INTEGER(GET_DIM(A))[0];
  long int nrhs=1;
  long int incx=1;
  long int incy=1;
  long int piv[n];
  
  SEXP tmpA, tmpb;
  
  PROTECT(tmpA = allocMatrix(REALSXP, n, n));
  PROTECT(tmpb = allocVector(REALSXP, n));
  long int nn = n*n;
  dcopy_(&nn, REAL(A), &incx, REAL(tmpA), &incy);
  dcopy_(&n, REAL(b), &incx, REAL(tmpb), &incy);
  
  printf("N = %ld\n", n);
  dgesv_(&n, &nrhs, REAL(tmpA), &n, piv, REAL(tmpb), &n, INTEGER(info));
  for(i=0; i<n; ++i) printf("%lf\n", REAL(b)[i]);
  UNPROTECT(2);
  return(tmpb);
}
  • RのAPIをCで使うときのお作法なのですが、新たにメモリを確保する場合はPROTECTというマクロを使わねばなりません.勝手にガベージコレクションされるのを防ぐためのようです.そして、それをどこかでUNPROTECTする必要があります.
  • 「dcopy_」はblasの関数で、変数の中身をコピーしてくれます.Aはn×nなので、二乗しています.
  • Rから呼ぶ
> A <- cbind(c(1, 3, 1), c(1, 1, -2), c(1, -3, -5))
> b <- c(3,1,-6)
> solve(A, b)
[1] 1 1 1
> info <- integer(1)
> .Call("Rdgesv01", A, b, info)
N = 3
3.000000   #元のbを表示している.
1.000000   #元のbを表示している.
-6.000000   #元のbを表示している.
[1] 1 1 1   # 結果はこの通り
> A   #書き換わっていない
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    3    1   -3
[3,]    1   -2   -5
> b   #書き換わっていない
[1]  3  1 -6
  • パチパチ

Cの中でRの関数を評価する

  • これを実現
print(x, digits=3)
  • C
SEXP eval_fun01(SEXP a, SEXP rho)
{
  SEXP s, t;
  SEXP ret;
  PROTECT(t = s = allocList(3));
  // Setting the type to LANGSXP makes this a call which can be evaluated
  SET_TYPEOF(s, LANGSXP);
  // set a symbol (pointing to the function to be called)
  SETCAR(t, install("print")); t = CDR(t);
  // set the first unnamed
  SETCAR(t, a); t = CDR(t);
  // set the second named
  SETCAR(t, allocVector(INTSXP, 1));
  INTEGER(CAR(t))[0] = 3;
  SET_TAG(t, install("digits"));
  PROTECT(ret = eval(s, rho));
  //PrintValue(ret);
  UNPROTECT(2);
  return R_NilValue;
}
  • Rで実行
> print("Hello", digits=3)
[1] "Hello"
> .Call("eval_fun01", "Hello", .GlobalEnv)
[1] "Hello"
NULL
> print(3.123456789, digits=3)
[1] 3.12
> .Call("eval_fun01", 3.123456789, .GlobalEnv)
[1] 3.12
NULL

Cの中でRの関数を評価する

  • 引数が1つだけなら「lang2」が使えます.
  • 「print("Hello")」を実現
SEXP eval_fun02(SEXP a, SEXP rho)
{
  SEXP call;
  SEXP ret;
  
  PROTECT(call = lang2(install("print"), a));
  eval(call, rho);
  UNPROTECT(1);
  return R_NilValue;
}
  • Rで実行
> .Call("eval_fun02", "Hello", .GlobalEnv)
[1] "Hello"
NULL

Cの中でRの関数を評価する

  • 評価したいexpressionを「quote」でわたす.
SEXP eval_fun03(SEXP fr, SEXP env)
{
  SEXP ans;
  
  if(!isEnvironment(env)) error(" 'env' should be an environment");
  PROTECT(ans = eval(fr, env));
  UNPROTECT(1);
  return(ans);
}
  • Rで実行
> f <- function(x) x^2
> .Call("eval_fun03", quote(f(x)), new.env())
 以下にエラー f(x) :  オブジェクト 'x' がありません
  • 「f(x)」を計算したいのですが、「new.env()」で指定した環境?にxがないので、xが見つからないといわれます.グローバル環境でxを作ってあげると、「new.env()」にxがない場合にはグローバル環境まで探しにきてくれます.
> f <- function(x) x^2
> .Call("eval_fun03", quote(f(x)), new.env())
 以下にエラー f(x) :  オブジェクト 'x' がありません
> x <- 2
> .Call("eval_fun03", quote(f(x)), new.env())
[1] 4
> x <- pi
> .Call("eval_fun03", quote(f(x)), new.env())
[1] 9.869604

Cの中でRの関数を評価する

  • をもうちょっと一般的にやろうとするとこんな感じになるでしょうか?
SEXP eval_fun04(SEXP fr, SEXP x, SEXP env)
{
  SEXP ans;
  
  if(!isEnvironment(env)) error(" 'env' should be an environment");
  defineVar(install("x"), x, env);
  PROTECT(ans = eval(fr, env));
  UNPROTECT(1);
  return(ans);
}
  • これはCの中で、渡された新環境でxを定義していますから、グローバル環境にxがなくてもエラーはでません.
> .Call("eval_fun03", quote(f(x)), new.env())
[1] 9.869604   #の例
> f <- function(x) x^3
> .Call("eval_fun04", quote(f(x)), 2, new.env())
[1] 8
  • でもx1だと・・・
> .Call("eval_fun04", quote(f(x1)), 2, new.env())
 以下にエラー f(x1) :  オブジェクト 'x1' がありません

Cの中でRの関数を評価する

  • 文字列をそのまま評価したいというそんなあなたには・・・
  • Cのソース
    • 【#include <R_ext/Parse.h>】の追加が必要です.
SEXP eval_fun05(SEXP a, SEXP env)
{
  SEXP cv, call, call2;
  SEXP ans;
  int i;
  ParseStatus status;
  char *cmd1 = "print(a)";  //なんでもOK.文字列を準備
  PROTECT(cv = allocVector(STRSXP, 1));  //文字列用のオブジェクトを用意
  SET_STRING_ELT(cv, 0, mkChar(cmd1));  //文字列変数に入れる
  PROTECT(call = R_ParseVector(cv, 1, &status, R_NilValue));  //お作法です.
  PROTECT(call2 = lang2(install("eval"), call));  //eval関数にセットする
  for(i=0; i<5; i++){
    INTEGER(a)[0] = INTEGER(a)[0] + i;
    defineVar(install("a"), a, env);
    PROTECT(ans = eval(call2, env));
    UNPROTECT(1);
  }
  UNPROTECT(3);
  return ans;
}
  • Rの結果
> dyn.load("test.so")
> is.loaded("eval_fun05")
[1] TRUE
> .Call("eval_fun05", as.integer(2), new.env())
[1] 2
[1] 3
[1] 5
[1] 8
[1] 12
  • 「call2」というオブジェクト(名前はなんでも良い)を用意しておけば、例えばループの中で何度も評価できる.
  • 評価する前に、(この例では)表示させる「a」というオブジェクトを「R_env」環境できちんと定義しておく必要がある.そうでなければ、何度評価しても、同じ値しか表示されません(笑
  • 文字列をRから渡すなり、Cにハードコーディングするなり、お好きにどうぞ.
  • 管理人みゅもこれを多様.

コメント