Posts Tagged ‘Dynamical System’
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
Rの練習のために力学系を描いてみた
メモをまとめてて思い出したんだけど
力学系使ったロボティクスって、他の分野からだと「エージェントの初期状態決めると後は決定論的に決まる」と受け止められてるらしい。
惜しい。ちょっと違うから注意して欲しい。
「系の初期状態決めると後は決定論的に決まる」だからね！
だから、エージェントの外の系が決定論的に（例えば微分方程式で）かけなかったり（例えば、時間による周期外力が加わったり）すると力学系ではない。もっと、環境も系の一部だって認識してもらえないんだろうか？
逆説的だけど、もっと環境に着目した力学系使ったお話あっても面白いと思った。
