![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
Rから他言語利用
(単なる雛型の提供です。誤りも含まれてると思います。修正お願いします。FortranやC++バージョンも拡充されるといいですね。)
ここでは、RからC, Fortranなど他言語で記述されたルーチン(関数)を利用する方法を説明します。 詳細は，R付属のドキュメント「Writing R Extensions」の「System and foreign language interfaces」の章に詳しいのですが，記述の順序が実際の作業手順と一致しておらず，分りにくいので，ここでは実際の作業手順に則して記述します。
例題として，単純なスカラーどうしの足し算をおこなう関数をCで書いて，そ れをRで利用することを試みる。
C 言語でルーチン(関数)を記述して，mysum.c というファイル名で保存する。 内容は以下のとおり。
void mysum(double *a, double *b, double *c)
{
double sum;
sum = *a + *b;
*c = sum;
}
ここで注意すべきは，
共用ライブラリをビルドするためにターミナルにて，以下のコマンドを実行す る。
$ R CMD SHLIB mysum.c
そうすると，オブジェクトファイル mysum.o 、共用ライブラリ mysum.so が作 成される。
R 内で共用ライブラリをロードするために以下のコマンドを実行する。
> dyn.load("./mysum.so")
ここを参考にして下さい。
必要なものは以下の三つで、ダウンロード後、展開してパスを通します。
C/Fortran のプログラム作成後（ここでは"mysum.c"）、コマンドプロンプトから
Rcmd SHLIB mysum.c
を実行、すると dll ファイルができるので、 R のコンソールから、
> dyn.load("mysum.dll")
として読み込みます。あとの手順は、以下と同様です。
Q&Aからなかまさんのありがたいお言葉を転載します。
上記mno-cygwin-howto.txtはリンク切れで見つからなくなっています。和訳がhttp://www.sixnine.net/cygwin/translation/devel/mno-cygwin-howto.htmlにあります。
R 内にて，関数 .C() を用いて関数 mysum() を定義する。
mysum <- function(a,b){
.C("mysum", #Cルーチン名
arg1=as.double(a), #Cルーチン第1引数の型指定
arg2=as.double(b), #Cルーチン第2引数
arg3=double(1) #Cルーチン第3引数
)
}
ここで注意すべきは，
(1).Cの第1引数には，C言語で記述したルーチン名(mysum)を文字列で指定
(2).Cの残りの引数は，ルーチンの引数と「順番」「型」「数」が対応するように列挙
(3)ここで定義した関数は，ルーチンの全引数を要素とするリストを返す。よって，
特定の引数のみ(例えば足し算の答)を返すようにしたければ，最後か ら2行目を
")$arg3"などとする。
上記の関数定義を行なったあとで，
> z <- mysum(1.234,2.345) > z
などとしてみると以下のようにな結果になる。
$arg1 [1] 1.234 $arg2 [1] 2.345 $arg3 [1] 3.579
パラメータベクトルを C サブルーティンに渡して、Mersenne-Twister 乱数を用いたモンテカルロ法を実行。結果の 100x200 整数行列をベクトルとして返す例 (実際に使ったものの抜粋)。
一回のシミュレーションを実行する R 関数
onesimulation <- function () {
cat("simulation start \n")
N <- GMRF(c(2.632, -0.9675, -0.2236, -0.3041, 5000))
dim(N) <- c(100,200) # ベクトル返り値を行列に変換
write(t(N), ncolumns=200, file="simulation.result")
}
与えられたパラメータベクトル param を用い、C サブルーティン gmrf を呼び出す R 関数
GMRF <- function (param) .Call( "gmrf", as.numeric(as.vector(param)) )
C サブルーティン gmrf.c は Mersenne-Twister 乱数発生 C サブルーティン mt19937-2.c とともにコンパイルしておく。
$ R CMD SHLIB gmrf.c mt19937-2.c
オブジェクトファイル gmrf.o 共用ライブラリgmrf.so が作成される。
R内で共用ライブラリをロードするために以下のコマンドを実行する。
> dyn.load("./gmrf.so")
/////////////////////////////////////////////////////////////
// シミュレーションプログラム gmrf.c
////////////////////////////////////////////////////////////
#include <R.h>
#include <Rinternals.h>
// use mt19937-2.c (mersenne-twister algorithm) Mersenne-Twister 乱数使用
extern void sgenrand(unsigned long seed);
extern double genrand(void);
//////////////////////////////////////////////////////////////////////////////////////
// 与えられたパラメータベクトル param でシミュレーション
//////////////////////////////////////////////////////////////////////////////////////
SEXP gmrf(SEXP param)
{
long int niter; // 繰り返し回数
double p[5]; // パラメータ用ベクトル
int N[100][200]; // 計算結果用整数行列
--- (その他省略) -----
SEXP ans; // 返り値用 R データの宣言
for(i = 0; i <= 3; i++) p[i] = REAL(param)[i]; // パラメータを取り出し
niter = REAL(param)[4]; // シミュレーション回数を取り出し
PROTECT(ans = allocVector(INTSXP, 100*200)); // 返り値用ベクトルを作成、保護
sgenrand(time(NULL)); // Mersenne-Twister 乱数初期化
--- (途中省略) -----
for(i = 0; i < 10; i++) {
(genrand() < 0.6) ? (Pat[i]=1):(Pat[i]=0); // Mersenne-Twister 乱数使用例
}
--- (途中省略) -----
for (i = 0; i < 100; i++)
for (j = 0; j < 200; j++)
INTEGER(ans)[200*i+j] = N[i][j]; // 結果をベクトル ans にコピー
UNPROTECT(1); // 変数 ans の PROTECT を解除
return(ans); // ベクトル返り値
}
Linux の場合、 zip ファイルを解凍し(unzip mandelbrot.zip)、コンソールで
$ R CMD SHLIB mandelbrot.c
とすると共有ライブラリ mandelbrot.so が出来る。R を起動し
> dyn.load("mandelbrot.so") # 共有ライブラリをロード
> source("mandelc.R") # R コードを読み込み
> image(mandelbrot()) # 実行と image 関数による表示(オプション引数で色々変化)
> mandelbrot <- # オプション引数で色々変化
function(x = c(-3, 1), # x limits
y = c(-1.8, 1.8), # y limits
nx = 600, # x resolutio
ny = 600, # y resolution
iter = 20) # maximun number of iterations
{関数本体省略}
Windowsの場合も、MinGWだと
$ Rcmd.exe SHLIB mandelbrot.c
で良いが、これを無理矢理Cygwinでコンパイルしたい場合は以下のようにする。
1. cc1.exe と libgcc.a の在り処を探す．私の環境では /usr/lib/gcc/i686-pc-cygwin/3.4.4 に双方ともあった（binutilsパッケージを選択していればインストールされるはず．他にもmingw関連のパッケージを選択しなければならない）
2. 上記ディレクトリにパスを通す
export PATH=$PATH:/usr/lib/gcc/i686-pc-cygwin/3.4.4
3. その上でソースをコンパイルする．-mno-cygwin オプションがミソ
gcc -shared -o mandelbrot.dll mandelbrot.c -mno-cygwin -L/usr/lib/gcc/i686-pc-cygwin/3.4.4 -lgcc
4. mandelbrot.dll が出来るので、dyn.loadすればO.K.
実際にはCygwinのバージョンにも依存しますし、なかまさんご指摘のように一筋縄ではいきません。MinGWを入れてRの枠組みでやる方が遥かに簡単です。
直角三角形の，直角を挟む二辺の長さを与え，斜辺の長さを返す関数の例。
function ではなく，subroutine にする。hyp.f という名前にした。拡張子は .f 戻り値は変数にセットする（例では変数 c）
subroutine hyp(a, b, c) real*8 a, b, c c = sqrt(a**2+b**2) end
hyp.so という名前のファイルができる
$ R CMD SHLIB hyp.f gfortran-4.0 -arch ppc -fPIC -fno-common -g -O2 -c hyp.f -o hyp.o gcc-4.0 -arch ppc -flat_namespace -bundle -undefined suppress -L/usr/local/lib -o hyp.so hyp.o -lgfortran -lgcc_s -lSystemStubs -lmx -lSystem -L/Library/Frameworks/R.framework/Resources/lib/ppc -lR
hyp <- function(a,b){
.Fortran("hyp",
arg1=as.double(a),
arg2=as.double(b),
arg3=double(1)
)
}
> dyn.load("hyp.so")
> hyp <- function(a,b){
+ .Fortran("hyp",
+ arg1=as.double(a),
+ arg2=as.double(b),
+ arg3=double(1)
+ )
+ }
> hyp(3,4)
$arg1
[1] 3
$arg2
[1] 4
$arg3
[1] 5
> hyp(1.2345, 3.1415)
$arg1
[1] 1.2345
$arg2
[1] 3.1415
$arg3
[1] 3.375354
