Rから他言語利用
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
//佐藤(2003-7-26)
COLOR(red){SIZE(30){Rから他言語利用}}
(単なる雛型の提供です。誤りも含まれてると思います。修正お願いします。FortranやC++バージョンも拡充されるといいですね。)
-CRAN にあるパッケージ RcppTemplate は R と Cpp 間のインタフェイスのテンプレートを提供する
-mno-cygwin-howto.txtの和訳へのリンクを追加したのと、マンデルブロ集合のCプログラムを無理矢理Cygwinでコンパイルするための段取りを追記しました。 -- &new{2005-12-21 (水) 16:19:39};
-FORTRAN を利用する場合を付け加えました（C とほとんど，同じですね） -- &new{2006-06-07 (水) 22:30:01};
- 私の環境(WinXP)でR2.5.0＆2.5.1のwin32でうまく行かずR2.3.1でできました。バグなのかよくわからないですが、初心者Q&Aに報告します。 -- [[たけ]] &new{2007-07-05 (木) 17:35:36};
- はやとちりだったようです。他の方はうまくいっていますので、最近のmingwがなにかよくなかったのかな？ -- [[たけ]] &new{2007-07-10 (火) 11:53:54};
- [[こんな:http://www.okada.jp.org/RWiki/?R%A4%AB%A4%E9%C2%BE%B8%C0%B8%EC%CD%F8%CD%D1%A1%CA%A5%A4%A5%F3%A5%C1%A5%AD%C8%C7%A1%CB]]インチキな方法があるのですね。大歓迎です -- &new{2007-08-07 (水) 20:11:21};
#comment
ここでは、RからC, Fortranなど他言語で記述されたルーチン(関数)を利用する方法を説明します。
詳細は，R付属のドキュメント「Writing R Extensions」の「System and foreign language interfaces」の章に詳しいのですが，記述の順序が実際の作業手順と一致しておらず，分りにくいので，ここでは実際の作業手順に則して記述します。
// "*", "**", "***" による段落に応じて、ハイパーリンク付きの目次一覧をこの位置に置く
// 行頭が "//" で始まる行はコメントで表示されません。
#contents
// SIZE(num){テキスト} はテキスト文字の大きさを指定します
// COLOR{color}{テキスト} はテキスト文字の色を指定します
// 行頭が "**" で始まると、直前の "*" で始まる段落の副段落として整形されます
// 二重括弧 "((...))" で囲ったテキストは注釈になります。注釈はページ末尾に置かれ、ハイパーリンクが張られます
// 行頭に半角空白を置くとこの通り枠で囲って見易く整形されます
// 行頭に半角空白がある引き続く行は一まとめに整形されます
// コードを書くのに便利です
// 行頭が "***" で始まると、直前の "**" で始まる副段落の副段落として整形されます// 行頭に "----" を置くとその場所に水平線が引かれます
----
// 行頭に "*" を置くと、段落見出しとして整形されます
*作業手順 [#wcb88f43]
+ C 言語でコード(ルーチン)を記述
+ C のコードから共用ライブラリをビルド
+ R 内にて共用ライブラリをロード
+ R 内で関数定義
+ 利用
*例題(C言語篇 その1) [#d5f86691]
例題として，単純なスカラーどうしの足し算をおこなう関数をCで書いて，そ
れをRで利用することを試みる。
** C 言語でコードを記述 [#g9f53826]
C 言語でルーチン(関数)を記述して，''mysum.c'' というファイル名で保存する。
内容は以下のとおり。
void mysum(double *a, double *b, double *c)
{
double sum;
sum = *a + *b;
*c = sum;
}
ここで注意すべきは，
+ ルーチン ''mysum'' は ''void型'' であること
+ 引数は，たとえスカラーであってもポインタとしておくこと (''double *'')
+ ほしい返り値（計算結果）も引数に含めておくこと (''double *c'')
** コードから共用ライブラリをビルド [#x715e8ca]
共用ライブラリをビルドするためにターミナルにて，以下のコマンドを実行す
る。
$ R CMD SHLIB mysum.c
そうすると，オブジェクトファイル ''mysum.o'' 、共用ライブラリ ''mysum.so'' が作
成される。
** R内で共用ライブラリをロード [#r5168c3c]
R 内で共用ライブラリをロードするために以下のコマンドを実行する。
> dyn.load("./mysum.so")
**Windows の場合 [#e81931ba]
[[ここ:http://gsbwww.uchicago.edu/fac/peter.rossi/teaching/37904/Adding%20Functions%20Written%20in%20C%20to%20R.pdf]]を参考にして下さい。
必要なものは以下の三つで、ダウンロード後、展開してパスを通します。~
+ A set of compatible tools which are various UNIX utilities ported to Windows~
http://www.murdoch-sutherland.com/Rtools/
+ A version of PERL for Windows~
http://www.activestate.com/Products/ActivePerl/Download.html.
+ The “minimal” version of GNU for windows (MinGW)~
http://prdownloads.sourceforge.net/mingw/~
(例えば，C:\MinGW\bin;C:\Perl\bin;C:\R\rw2001\binを環境変数のPATHに追加．tools.zipの中身は..\R\rw2001\bin下に置く)
C/Fortran のプログラム作成後（ここでは"mysum.c"）、コマンドプロンプトから
Rcmd SHLIB mysum.c
を実行、すると dll ファイルができるので、 R のコンソールから、
> dyn.load("mysum.dll")
として読み込みます。あとの手順は、以下と同様です。
Q&Aからなかまさんのありがたいお言葉を転載します。
-目的がRと他言語(C,Fortran)ならMingwを使いましょう.目的がCygwin上でなら http://www.nanotech.wisc.edu/~khan/software/gnu-win32/mno-cygwin-howto.txt が役立つ化もしれません. ただ,ヘッダーファイルの差異や面倒な事に悩んだり, 後でパッケージ化したり等不便を強いられると思います. R CMD SHLIB(Rの枠組を使うなら) Mingwを入れるべきですし, Cygwinを使うならR CMD SHLIBを使うべきでは無いと思います. 幸いな事に外部プログラムの作成に必要な環境のドキュメントは存在します.Rの外部プログラム作成に必要な環境とは、即ちRのビルド環境を揃える事です. -- [[なかま]] {2005-05-01 (日) 00:46:33};
上記mno-cygwin-howto.txtはリンク切れで見つからなくなっています。和訳がhttp://www.sixnine.net/cygwin/translation/devel/mno-cygwin-howto.htmlにあります。
** R 内で関数定義 [#o03829b7]
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"などとする。
**利用 [#hcc1e8bc]
上記の関数定義を行なったあとで，
> z <- mysum(1.234,2.345)
> z
などとしてみると以下のようにな結果になる。
$arg1
[1] 1.234
$arg2
[1] 2.345
$arg3
[1] 3.579
*例題(C言語篇 その 2) もう少し大規模な例 [#y28f32ad]
パラメータベクトルを 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); // ベクトル返り値
}
* Mandelbrot フラクタル図形を書くRプログラム例 (高速化のためCサブルーティンを利用している) 2003.12.15 [#s9f649e3]
[[ソースと実行例はここ:http://www.geocities.com/mariodosreis/Rfractals/]]
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の枠組みでやる方が遥かに簡単です。
*FORTRAN の場合（Mac OS X） Jun 07, 2006 [#eb80082a]
直角三角形の，直角を挟む二辺の長さを与え，斜辺の長さを返す関数の例。
** Fortran プログラムを書く [#t8e11e66]
function ではなく，subroutine にする。hyp.f という名前にした。拡張子は .f
戻り値は変数にセットする（例では変数 c）
subroutine hyp(a, b, c)
real*8 a, b, c
c = sqrt(a**2+b**2)
end
** 共用ライブラリを作る [#d50bf7ca]
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
** R プログラムを書く [#x5a447bd]
hyp <- function(a,b){
.Fortran("hyp",
arg1=as.double(a),
arg2=as.double(b),
arg3=double(1)
)
}
** R コンソールで利用する [#m92fa346]
> 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
*参考文献 [#y8641476]
+R Development Core Team,"Writing R Extensions",chapter 4.(R付属マニュアル)
注意点:実際の作業手順と記述順序が逆転しているので見通しが悪く分りにくい。(と思う)
+R.A.ベッカー+J.M.チェンバース+A.R.ウィルクス 著(渋谷政昭+柴田里程 訳)「S言語 データ解析とグラフィックスのためのプログラミング環境 I」共立出版，1991年，第7章。
注意点:この書は，S(engine)用の記述なので，手順自体は同じだが，R(engine)とは細かい点では異なる。例えば，Rではオブジェクトファイル(***.o)ではなく共用ライブラリ(***.so)をロードする点など。
+W.N.Venables and B.D.Ripley(2000),"S Programming",Springer-Verlag,chapter6，およびAppendix 2.
関連リンク(R,Cのコード) http://www.stats.ox.ac.uk/pub/MASS3/Sprog/index.html
+J.M. Chambers 著(垂水ら訳)「データによるプログラミング」森北出版，第11章。
+[[Rのお部屋（あーるのおへや）>http://phoenixx.sakura.ne.jp/R/wiki.cgi?page=FrontPage]]
--[[RからCを呼ぶ（Linux編）>http://phoenixx.sakura.ne.jp/R/wiki.cgi?page=BugTrack%2D%A3%D2%C8%F7%CB%BA%CF%BF%2F56]]
--[[RからCを呼ぶ（Linux編） - ちょっとしたコツ>http://phoenixx.sakura.ne.jp/R/wiki.cgi?page=BugTrack%2D%A3%D2%C8%F7%CB%BA%CF%BF%2F57]]
--[[RからCを呼ぶ（Linux編） - blas・lapack>http://phoenixx.sakura.ne.jp/R/wiki.cgi?page=BugTrack%2D%A3%D2%C8%F7%CB%BA%CF%BF%2F58]]
終了行:
//佐藤(2003-7-26)
COLOR(red){SIZE(30){Rから他言語利用}}
(単なる雛型の提供です。誤りも含まれてると思います。修正お願いします。FortranやC++バージョンも拡充されるといいですね。)
-CRAN にあるパッケージ RcppTemplate は R と Cpp 間のインタフェイスのテンプレートを提供する
-mno-cygwin-howto.txtの和訳へのリンクを追加したのと、マンデルブロ集合のCプログラムを無理矢理Cygwinでコンパイルするための段取りを追記しました。 -- &new{2005-12-21 (水) 16:19:39};
-FORTRAN を利用する場合を付け加えました（C とほとんど，同じですね） -- &new{2006-06-07 (水) 22:30:01};
- 私の環境(WinXP)でR2.5.0＆2.5.1のwin32でうまく行かずR2.3.1でできました。バグなのかよくわからないですが、初心者Q&Aに報告します。 -- [[たけ]] &new{2007-07-05 (木) 17:35:36};
- はやとちりだったようです。他の方はうまくいっていますので、最近のmingwがなにかよくなかったのかな？ -- [[たけ]] &new{2007-07-10 (火) 11:53:54};
- [[こんな:http://www.okada.jp.org/RWiki/?R%A4%AB%A4%E9%C2%BE%B8%C0%B8%EC%CD%F8%CD%D1%A1%CA%A5%A4%A5%F3%A5%C1%A5%AD%C8%C7%A1%CB]]インチキな方法があるのですね。大歓迎です -- &new{2007-08-07 (水) 20:11:21};
#comment
ここでは、RからC, Fortranなど他言語で記述されたルーチン(関数)を利用する方法を説明します。
詳細は，R付属のドキュメント「Writing R Extensions」の「System and foreign language interfaces」の章に詳しいのですが，記述の順序が実際の作業手順と一致しておらず，分りにくいので，ここでは実際の作業手順に則して記述します。
// "*", "**", "***" による段落に応じて、ハイパーリンク付きの目次一覧をこの位置に置く
// 行頭が "//" で始まる行はコメントで表示されません。
#contents
// SIZE(num){テキスト} はテキスト文字の大きさを指定します
// COLOR{color}{テキスト} はテキスト文字の色を指定します
// 行頭が "**" で始まると、直前の "*" で始まる段落の副段落として整形されます
// 二重括弧 "((...))" で囲ったテキストは注釈になります。注釈はページ末尾に置かれ、ハイパーリンクが張られます
// 行頭に半角空白を置くとこの通り枠で囲って見易く整形されます
// 行頭に半角空白がある引き続く行は一まとめに整形されます
// コードを書くのに便利です
// 行頭が "***" で始まると、直前の "**" で始まる副段落の副段落として整形されます// 行頭に "----" を置くとその場所に水平線が引かれます
----
// 行頭に "*" を置くと、段落見出しとして整形されます
*作業手順 [#wcb88f43]
+ C 言語でコード(ルーチン)を記述
+ C のコードから共用ライブラリをビルド
+ R 内にて共用ライブラリをロード
+ R 内で関数定義
+ 利用
*例題(C言語篇 その1) [#d5f86691]
例題として，単純なスカラーどうしの足し算をおこなう関数をCで書いて，そ
れをRで利用することを試みる。
** C 言語でコードを記述 [#g9f53826]
C 言語でルーチン(関数)を記述して，''mysum.c'' というファイル名で保存する。
内容は以下のとおり。
void mysum(double *a, double *b, double *c)
{
double sum;
sum = *a + *b;
*c = sum;
}
ここで注意すべきは，
+ ルーチン ''mysum'' は ''void型'' であること
+ 引数は，たとえスカラーであってもポインタとしておくこと (''double *'')
+ ほしい返り値（計算結果）も引数に含めておくこと (''double *c'')
** コードから共用ライブラリをビルド [#x715e8ca]
共用ライブラリをビルドするためにターミナルにて，以下のコマンドを実行す
る。
$ R CMD SHLIB mysum.c
そうすると，オブジェクトファイル ''mysum.o'' 、共用ライブラリ ''mysum.so'' が作
成される。
** R内で共用ライブラリをロード [#r5168c3c]
R 内で共用ライブラリをロードするために以下のコマンドを実行する。
> dyn.load("./mysum.so")
**Windows の場合 [#e81931ba]
[[ここ:http://gsbwww.uchicago.edu/fac/peter.rossi/teaching/37904/Adding%20Functions%20Written%20in%20C%20to%20R.pdf]]を参考にして下さい。
必要なものは以下の三つで、ダウンロード後、展開してパスを通します。~
+ A set of compatible tools which are various UNIX utilities ported to Windows~
http://www.murdoch-sutherland.com/Rtools/
+ A version of PERL for Windows~
http://www.activestate.com/Products/ActivePerl/Download.html.
+ The “minimal” version of GNU for windows (MinGW)~
http://prdownloads.sourceforge.net/mingw/~
(例えば，C:\MinGW\bin;C:\Perl\bin;C:\R\rw2001\binを環境変数のPATHに追加．tools.zipの中身は..\R\rw2001\bin下に置く)
C/Fortran のプログラム作成後（ここでは"mysum.c"）、コマンドプロンプトから
Rcmd SHLIB mysum.c
を実行、すると dll ファイルができるので、 R のコンソールから、
> dyn.load("mysum.dll")
として読み込みます。あとの手順は、以下と同様です。
Q&Aからなかまさんのありがたいお言葉を転載します。
-目的がRと他言語(C,Fortran)ならMingwを使いましょう.目的がCygwin上でなら http://www.nanotech.wisc.edu/~khan/software/gnu-win32/mno-cygwin-howto.txt が役立つ化もしれません. ただ,ヘッダーファイルの差異や面倒な事に悩んだり, 後でパッケージ化したり等不便を強いられると思います. R CMD SHLIB(Rの枠組を使うなら) Mingwを入れるべきですし, Cygwinを使うならR CMD SHLIBを使うべきでは無いと思います. 幸いな事に外部プログラムの作成に必要な環境のドキュメントは存在します.Rの外部プログラム作成に必要な環境とは、即ちRのビルド環境を揃える事です. -- [[なかま]] {2005-05-01 (日) 00:46:33};
上記mno-cygwin-howto.txtはリンク切れで見つからなくなっています。和訳がhttp://www.sixnine.net/cygwin/translation/devel/mno-cygwin-howto.htmlにあります。
** R 内で関数定義 [#o03829b7]
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"などとする。
**利用 [#hcc1e8bc]
上記の関数定義を行なったあとで，
> z <- mysum(1.234,2.345)
> z
などとしてみると以下のようにな結果になる。
$arg1
[1] 1.234
$arg2
[1] 2.345
$arg3
[1] 3.579
*例題(C言語篇 その 2) もう少し大規模な例 [#y28f32ad]
パラメータベクトルを 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); // ベクトル返り値
}
* Mandelbrot フラクタル図形を書くRプログラム例 (高速化のためCサブルーティンを利用している) 2003.12.15 [#s9f649e3]
[[ソースと実行例はここ:http://www.geocities.com/mariodosreis/Rfractals/]]
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の枠組みでやる方が遥かに簡単です。
*FORTRAN の場合（Mac OS X） Jun 07, 2006 [#eb80082a]
直角三角形の，直角を挟む二辺の長さを与え，斜辺の長さを返す関数の例。
** Fortran プログラムを書く [#t8e11e66]
function ではなく，subroutine にする。hyp.f という名前にした。拡張子は .f
戻り値は変数にセットする（例では変数 c）
subroutine hyp(a, b, c)
real*8 a, b, c
c = sqrt(a**2+b**2)
end
** 共用ライブラリを作る [#d50bf7ca]
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
** R プログラムを書く [#x5a447bd]
hyp <- function(a,b){
.Fortran("hyp",
arg1=as.double(a),
arg2=as.double(b),
arg3=double(1)
)
}
** R コンソールで利用する [#m92fa346]
> 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
*参考文献 [#y8641476]
+R Development Core Team,"Writing R Extensions",chapter 4.(R付属マニュアル)
注意点:実際の作業手順と記述順序が逆転しているので見通しが悪く分りにくい。(と思う)
+R.A.ベッカー+J.M.チェンバース+A.R.ウィルクス 著(渋谷政昭+柴田里程 訳)「S言語 データ解析とグラフィックスのためのプログラミング環境 I」共立出版，1991年，第7章。
注意点:この書は，S(engine)用の記述なので，手順自体は同じだが，R(engine)とは細かい点では異なる。例えば，Rではオブジェクトファイル(***.o)ではなく共用ライブラリ(***.so)をロードする点など。
+W.N.Venables and B.D.Ripley(2000),"S Programming",Springer-Verlag,chapter6，およびAppendix 2.
関連リンク(R,Cのコード) http://www.stats.ox.ac.uk/pub/MASS3/Sprog/index.html
+J.M. Chambers 著(垂水ら訳)「データによるプログラミング」森北出版，第11章。
+[[Rのお部屋（あーるのおへや）>http://phoenixx.sakura.ne.jp/R/wiki.cgi?page=FrontPage]]
--[[RからCを呼ぶ（Linux編）>http://phoenixx.sakura.ne.jp/R/wiki.cgi?page=BugTrack%2D%A3%D2%C8%F7%CB%BA%CF%BF%2F56]]
--[[RからCを呼ぶ（Linux編） - ちょっとしたコツ>http://phoenixx.sakura.ne.jp/R/wiki.cgi?page=BugTrack%2D%A3%D2%C8%F7%CB%BA%CF%BF%2F57]]
--[[RからCを呼ぶ（Linux編） - blas・lapack>http://phoenixx.sakura.ne.jp/R/wiki.cgi?page=BugTrack%2D%A3%D2%C8%F7%CB%BA%CF%BF%2F58]]
ページ名:
![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
