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


!!!R（R言語）からCを呼ぶ（Linux編） - ちょっとしたコツ
*投稿者： みゅ
*カテゴリ： なし
*優先度： 普通
*状態： 完了
*日時： 2010年01月23日 20時41分35秒
//{{bugstate}}
!!内容
*「ちょっとしたコツ」集
*管理人みゅの備忘録（笑
*ここから先は、[[BugTrack-Ｒ備忘録/56]]のCのソースにどんどん追加していくだけで、動くはずです．
*ソース → {{ref test.c}}
*バグ報告や質問は[別館|http://phoenixx.sakura.ne.jp/mt/R/2010/01/rc.html]でお願いします．
!!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から、Aｘ=ｂを解く
 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の書き方だと、引数にした「ｂ」自体の中身も書き換わってしまっています．dgesvの仕様なので仕方ありません．（ついでに言うとAも書き換わる・・・）
 > b
 [1] 1 1 1
*これを、防ぐには結果を入れる場所を用意して、そこに結果をセットして返すしかありません．
*Rでこれを実現するのは簡単です．「A」「ｂ」が既にセットされていると仮定すると
 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はｎ×ｎなので、二乗しています．
*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   #元のｂを表示している．
 1.000000   #元のｂを表示している．
 -6.000000   #元のｂを表示している．
 [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の関数を評価する②
*引数が１つだけなら「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()」で指定した環境？にｘがないので、ｘが見つからないといわれます．グローバル環境でｘを作ってあげると、「new.env()」にｘがない場合にはグローバル環境まで探しにきてくれます．
 > 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の中で、渡された新環境でｘを定義していますから、グローバル環境にｘがなくてもエラーはでません．
 > .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にハードコーディングするなり、お好きにどうぞ．
*管理人みゅもこれを多様．


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