ナンクル力学系

学んだ事を書き連ねていこう。

Rで離散時間の力学系を計算するためにC拡張を書いた

leave a comment »

x_0 を与えたときに、 x_{t+1} = f(x_t) を使って 時系列 (x_0, x_1, x_2, \dots ) を計算するって単純なことが したいだけだったのに、上手くRでやる方法が思いつかなかった。

ので、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)&#91;time*i&#93; = REAL(x0)&#91;i&#93;;
  }
  for(t = 1; t < time; t++) {
    PROTECT( x1 = allocVector(REALSXP, dim) );
    for (i = 0; i < dim; ++i){
      REAL(x1)&#91;i&#93; = REAL(xt)&#91;t-1 + time*i&#93;;
    }
    defineVar(install("x"), x1, rho);
    UNPROTECT(1);

    PROTECT( x2 = coerceVector(eval(fn, rho), REALSXP) );
    for(i = 0; i < dim; i++){
      REAL(xt)&#91;t + time*i&#93; = REAL(x2)&#91;i&#93;;
    }
    UNPROTECT(1);
  }
  UNPROTECT(1);
  return(xt);
}
&#91;/sourcecode&#93;</div>
<div class="outline-5">
<h5>Rの関数(ラッパ?):</h5>

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) )
}
&#91;/sourcecode&#93;</div>
<div class="outline-5">
<h5>実行例:</h5>

> 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&#91;t-1,&#93;) })
   ユーザ   システム       経過
     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)&#91;time*i&#93; = REAL(x0)&#91;i&#93;;
  }
  for(t = 1; t < time; t++) {
    for (i = 0; i < dim; ++i){
      REAL(x1)&#91;i&#93; = REAL(xt)&#91;t-1 + time*i&#93;;
    }

    PROTECT( x2 = coerceVector(eval(fn, rho), REALSXP) );
    for(i = 0; i < dim; i++){
      REAL(xt)&#91;t + time*i&#93; = REAL(x2)&#91;i&#93;;
    }
    UNPROTECT(1);
  }
  UNPROTECT(2);
  return(xt);
}
&#91;/sourcecode&#93;</div>
<div class="outline-5">
<h5>実行例:</h5>

> 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&#91;t-1,&#93;) })
   ユーザ   システム       経過
     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
Advertisements

Written by tkf

March 11, 2009 at 7:17 pm

Posted in PC

Tagged with , ,

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: