Posts Tagged ‘C++’
Rで離散時間の力学系を計算するためにC拡張を書いた
ので、Cで拡張作った。
とりあえず動いたver
Cのコード
#include <R.h>
#include <Rinternals.h>
SEXP dtds(SEXP _time, SEXP x0, SEXP fn, SEXP rho)
{
int i, dim = length(x0);
int t, time = INTEGER(_time)[0];
SEXP xt, x1, x2;
PROTECT(xt = allocMatrix(REALSXP, time, dim));
for(i = 0; i < dim; i++){
REAL(xt)[time*i] = REAL(x0)[i];
}
for(t = 1; t < time; t++) {
PROTECT( x1 = allocVector(REALSXP, dim) );
for (i = 0; i < dim; ++i){
REAL(x1)[i] = REAL(xt)[t-1 + time*i];
}
defineVar(install("x"), x1, rho);
UNPROTECT(1);
PROTECT( x2 = coerceVector(eval(fn, rho), REALSXP) );
for(i = 0; i < dim; i++){
REAL(xt)[t + time*i] = REAL(x2)[i];
}
UNPROTECT(1);
}
UNPROTECT(1);
return(xt);
}
Rの関数（ラッパ？）:
dyn.load("dtds.so")
dtds <- function(N, x0, fmap){
.Call("dtds", as.integer(N), as.double(x0), quote(fmap(x)), new.env())
}
fmap <- function(x){
x * (x - 1)
}
rms <- function(x,y){
d <- c(x)-c(y)
sqrt( d %*% d/length(d) )
}
実行例:
> N <- 100000
> x0 <- c(-0.1,0.1)
>
> xt <- matrix(rep(x0, N),N,byrow=T)
> system.time(for (t in 2:N){ xt[t,] <- fmap(xt[t-1,]) })
ユーザ システム 経過
1.936 0.040 1.977
> print(xt[N,])
[1] 0.002237959 -0.002232851
>
> system.time(xs <- dtds(N, x0, fmap))
ユーザ システム 経過
0.42 0.00 0.42
# for使って計算したやつの4倍くらい?
> print(xs[N,])
[1] 0.002237959 -0.002232851
>
> print(rms(xt,xs))
[,1]
[1,] 0
# for使って計算したやつと同じ結果！
あれ？
4倍しかはやくならない？
くやしかったので少し改良したver
- 毎ステップ allocVector しなくて良いはず
- defineVar もボトルネックか？
Cのコード
#include <R.h>
#include <Rinternals.h>
SEXP dtds(SEXP _time, SEXP x0, SEXP fn, SEXP rho)
{
int i, dim = length(x0);
int t, time = INTEGER(_time)[0];
SEXP xt, x1, x2;
PROTECT(xt = allocMatrix(REALSXP, time, dim));
PROTECT( x1 = allocVector(REALSXP, dim) );
defineVar(install("x"), x1, rho);
for(i = 0; i < dim; i++){
REAL(xt)[time*i] = REAL(x0)[i];
}
for(t = 1; t < time; t++) {
for (i = 0; i < dim; ++i){
REAL(x1)[i] = REAL(xt)[t-1 + time*i];
}
PROTECT( x2 = coerceVector(eval(fn, rho), REALSXP) );
for(i = 0; i < dim; i++){
REAL(xt)[t + time*i] = REAL(x2)[i];
}
UNPROTECT(1);
}
UNPROTECT(2);
return(xt);
}
実行例:
> N <- 100000
> x0 <- c(-0.1,0.1)
> #x0 <- c(0.1)
>
> xt <- matrix(rep(x0, N),N,byrow=T)
> system.time(for (t in 2:N){ xt[t,] <- fmap(xt[t-1,]) })
ユーザ システム 経過
1.944 0.040 1.981
> #print(length(xt[1,]))
> #print(length(xt[,1]))
> print(xt[N,])
[1] 0.002237959 -0.002232851
>
> system.time(xs <- dtds(N, x0, fmap))
ユーザ システム 経過
0.388 0.004 0.392
> #print(length(xs[1,]))
> #print(length(xs[,1]))
> print(xs[N,])
[1] 0.002237959 -0.002232851
>
> print(rms(xt,xs))
[,1]
[1,] 0
うーん
- ほとんど変化なし
- やはり、evalがだめな感じがする
感想
PythonのC拡張比べてみて。（昔書いたPythonのC拡張はこんな感じ。 この例と同じ離散時間の系の計算だけど、fmapにあたる写像は決め打ち。）
- PythonでC拡張書くより、全然書く量が少ないのは良い
- ドキュメントがなさすぎるw
- Writing R Extensions だけ
- リファレンスが無いので、関数の役目が謎なままつかうことに
- 関数とかマクロとかフリーダム
- PythonのCコードは必ず「Py」から始まっていたり、マクロならそれ以下が大文字 になっていたりと安心する
- まー、Cだし、コンパイラが怒るので問題なし？
- 何回でも dyn.load(“dtds.so”) と dyn.unload(“dtds.so”) できるのはいいね！
- Pythonだと、一度読み込んだCのモジュールは reload しても変化なしなので，Python自体を再起動する必要があった
- デバッグ楽
- Cコードレベルのミス（？）に気づくRいいね
- 何度、 「 エラー： VECTOR_ELT() can only be applied to a ‘list’, not a ‘double’」 と怒られたことか。。。
- 慣れたらサクサクかけそう！
- 慣れる ＝ Rのソース読めるくらいか？w
PythonのC拡張でAPIを公開する例
PythonのC拡張を普通に書くと、その関数などは他のC拡張モジュールから使えない。という訳で他のC拡張から読み込めるようにAPIを公開する必要がある訳だけど、その作り方が少しトリッキーだったので公開してみる。
参考：拡張モジュールに C API を提供すると、NumPyのソース。
（追記。そういえばPython/C拡張の練習で書いたコードがあったので晒す：http://snipplr.com/view/8215/example-of-numpyc-api/ C APIは使ってないけど。）
あんまりPythonのC拡張どころか、C言語もそんなに触って無いのでおかしい所があるかも。ツッコミ大歓迎です！
準備するのは以下のファイル。son.Sonとfather.Fatherというモジュールを作ってFatherのメンバにSonをもってきてみた。
STLのシャッフルアルゴリズムstd::random_shuffleの乱数にBoost/random.hppを使う
C++の標準ライブラリSTLに入ってるシャッフルを使いたいんだけど，ちゃんとした乱数を使いたいってことでboostを使ってみることにした．
コードはいつものようにsnipplrに：
http://snipplr.com/view/5907/stdrandomshuffle-and-boostrandomhpp/
C++はちゃんと勉強せずにプログラムを書いたり読んだりしてるだけなので，変な所があるかもしれない．
参考（というか↓を組み合わせたら↑になるんだけどねｗ）：
