ナンクル力学系

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

Archive for March 2009

今の時期に読む本は自分の研究の個性を決めることになりそうだと思いながら買った

leave a comment »

洋物ばっかりである.

Thermodynamics of Chaotic Systems: An Introduction

数学的形式としての統計力学の話を学ぶための本(候補) で買うぞと言って買ったので,届いてから一ヶ月くらいたっているかも. しかし,3月は忙しかった(←言い訳か!くそっ)のであまり進んでいない.1/4くらいか. 今は熱力学的な話の部分で,統計力学の知識も熱力の知識もないのに進める.けど,ちょっと 遅い感じは否めないか.力学系と情報理論の部分はまあまあすんなりいけた(はず). 一度,統計力学の教科書の初めの部分読むべきか.でも読むのが遅いのは3月忙しかt(ry

Introduction to the Theory of Statistics

2008年にお勧めだった本 – Seeking for my unique color. で syou先生が 薦めていた本.何回か買おうかと迷ったけど,とうとう買ってしまった.どうしよう.

Principles of Neural Science

The Swingy Brain: 神経回路のサブネットワーク  のコメント で神経科学やるなら,と薦められ,買おうか何回も迷っていた本. カンデルでググると日本語の資料も沢山ある. 通称カンデルと呼ばれるめちゃくちゃ有名な本らしい.

博学的な知識を要求するような学問は好きじゃないんだけど,脳に関係する研究するなら やっぱり基礎的な知識が無いと論文さえ読めないのでどうにかしなくちゃいけない. 日本語の名前覚えても意味なさそうなので,ラボにある神経科学ではなくてこっちを読もう.

The Nonlinear Workbook

Chaos, Fractals, Cellular Automata, Neural Networks, Genetic Algorithms, Gene Expression Programming, Support Vector Machine, Wavelets, Hidden Markov Models, Fuzzy Logic with C++, Java and SymbolicC++ Programs

が本の中身らしい.解説 + サンプルプログラムという形で,ここまで詰め込んだ本はそうそう 無いだろう,と思って買った.新宿の本屋をぶらぶらしていたら見つけて気になってた本. amazon.comでの評判はそこまで良くなさそうだった(解説に関して)けど,なんちゃってな解説 でも良いからまとまっているのはうれしい,はず. 2年前に手に入れたかったなあ!

Written by tkf

March 30, 2009 at 10:33 pm

Posted in 物理, 数学

Tagged with ,

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

Written by tkf

March 11, 2009 at 7:17 pm

Posted in PC

Tagged with , ,

Rで複数の3D曲線を描く関数を少し直した

leave a comment »

Rで複数の3D曲線を描く関数を作ってみたけど異様に遅い « ナンクル力学系 で作った関数を少し変更。

関数と実行例:

lines3d <- function(linelist, theta=25, phi=30, expand=1, xlab="X", ylab="Y", zlab="Z", ticktype="detailed", xlim=NULL, ylim=NULL, zlim=NULL ){ # x,y,z の上限と下限を求める limit <- matrix(c(range(linelist[[1]][,1]) , range(linelist[[1]][,2]) , range(linelist[[1]][,3]) ), c(3,2), byrow=T) for (lin in linelist) { for (i in seq(3)) { limit[i,] <- range( lin[,i], limit[i,] ) } } # limitが指定されていれば上書き if ( !is.null(xlim) ){ limit[1,] <- xlim } if ( !is.null(ylim) ){ limit[2,] <- ylim } if ( !is.null(zlim) ){ limit[3,] <- zlim } # 上限と下限の差が非零かチェック for (i in seq(3)) { if ( limit[i,1] == limit[i,2] ){ limit[i,1] <- limit[i,1] - 1 limit[i,2] <- limit[i,2] + 1 } } pmat <- persp(0:1, 0:1, matrix(,2,2), xlim=limit[1,], ylim=limit[2,], zlim=limit[3,], theta=theta, phi=phi, expand=expand, xlab=xlab, ylab=ylab, zlab=zlab, ticktype=ticktype ) for (lin in linelist) { lines(trans3d(lin[,1], lin[,2], lin[,3], pmat)) } return(pmat) } x0 <- seq(-1,1,length=1000) y0 <- cos(x0*6*pi) z0 <- rep(0,length(x0)) x1 <- seq(-1,1,length=200) y1 <- rep(0,length(x1)) z1 <- sin(x1*2*pi) lines3d(list(cbind(x0,y0,z0), cbind(x1,y1,z1), cbind(x1,y1+1,z1*0.2), cbind(x1,y1-1,z1*0.5))) [/sourcecode]

「異様に遅い」は嘘でした><

id:syou6162さんに system.time を教えてもらって 関数を実行する部分だけ時間はかってみたら、全然遅くなかった。

疑問

  • 「system.time(lines3d <- function(略))」 とやってみた
    • 実時間でかなり時間かかる
    • でも、 system.time が表示する実行時間は0秒
    • なぜ? Rは関数の評価に時間かるんだろうか?
  • 「source(‘ファイルへのパス’)」で実行してみると瞬殺
    • RじゃなくてESSの問題とか?

他やること

  • まじめにやるなら、
    • lineごとに色や線種を変える
    • 凡例を出す
    • 外で定義したpmatを使えるようにする
    • 。。。とかやるべきだけど、[あとでやる]
      そんなに複雑なグラフは今は書かないはず。

Written by tkf

March 11, 2009 at 2:15 am

Posted in PC

Tagged with

Rで複数の3D曲線を描く関数を作ってみたけど異様に遅い

leave a comment »

球の三次元プロット (r-help 記事より, 2003.12.27) << グラフィックス参考実例集:自作グラフィックス投稿欄 – RjpWiki を参考にして、作ってみた。

関数と実行例:

# # デフォルトで定義されているっぽい?
# trans3d <- function(x,y,z, pmat) { # tr <- cbind(x,y,z,1) %*% pmat # list(x = tr[,1]/tr[,4], y= tr[,2]/tr[,4]) # } lines3d <- function(linelist, theta=25, phi=30, expand=1, xlab="X", ylab="Y", zlab="Z", ticktype="detailed", xlim=NULL, ylim=NULL, zlim=NULL ){ # x,y,z の上限と下限を求める limit <- list(c( min(linelist[[1]][[1]]), max(linelist[[1]][[1]]) ) , c( min(linelist[[1]][[2]]), max(linelist[[1]][[2]]) ) , c( min(linelist[[1]][[3]]), max(linelist[[1]][[3]]) ) ) for (lin in linelist) { for (i in seq(length(limit))) { limit[[i]] <- c(min(min(lin[[i]]), limit[[i]][1]) , max(max(lin[[i]]), limit[[i]][2]) ) } } # # 上書き(この辺に,xlim=NULLあたりの処理をかくつもり) # 上限と下限の差が0より大きいかチェック for (i in seq(length(limit))) { if ( limit[[i]][1] == limit[[i]][2] ){ limit[[i]][1] <- limit[[i]][1] - 1 limit[[i]][2] <- limit[[i]][2] + 1 } } pmat <- persp(0:1, 0:1, matrix(,2,2), xlim=limit[[1]], ylim=limit[[2]], zlim=limit[[3]], theta=theta, phi=phi, expand=expand, xlab=xlab, ylab=ylab, zlab=zlab, ticktype=ticktype ) for (lin in linelist) { lines(trans3d(lin[[1]], lin[[2]], lin[[3]], pmat)) } return(pmat) } x0 <- seq(-1,1,length=1000) y0 <- cos(x0*6*pi) z0 <- rep(0,length(x0)) x1 <- seq(-1,1,length=200) y1 <- rep(0,length(x1)) z1 <- sin(x1*2*pi) lines3d(list(list(x0,y0,z0), list(x1,y1,z1))) [/sourcecode]

なんで遅い?

  • たぶん、
    • forを使いまくっている
      • Rっぽくないよな
    • listを使いまくっているデータ構造が、そもそもRっぽい処理に向かない?
    • linelistの一つ目の要素を二回しらべている
      • [あとでなおす] 手を抜きました><

Written by tkf

March 10, 2009 at 11:52 pm

Posted in PC

Tagged with

Tsukuba.R#4のLTに使ったRのソースとかやっとUPした。あと+何か。

leave a comment »

発表スライドは Tsukuba.R#4で力学系の分岐現象についてLTしてきた « ナンクル力学系 にある。

なぜかgithubにまとめておくとにした

+何か

今更なんだよ、と思うけれどw

  • Tsukuba.RのWikiにトラバ飛ばすのかなと思ったら俺だけだったので悲しかった。
    • id:syou6162さんがまとめたりしてるけど、もれがないように自主的に送るのが 良いのではと思ったりする
    • なんで俺は同じ記事でふたつも飛ばしてるんだ
  • 聞きたい人がいるような内容ではなかった(気がする)ので、次やるとしたらニーズが ある内容をやりたい
  • ふらふらしている時に良い発表はできません
    • 元気がなかったのでスライド飛ばした(ぉ
  • Tsukuba.Rの運営に関してもらった反応に応えていくよ!! – Seeking for my unique color.
    とそのリンク先を読んで、

    • 「Rというツールを通して他の分野のこと知る」みたいな感じになれば面白いなあ
    • まったく違う分野でも数学的な構造が同じ、ってことは良くある。この「数学的な構造」を 「Rの使い方」に置き換えても成り立つはず
    • これは、Rの使い方(テクニック)を披露しあう場になれば当然成り立つはず
    • だから「背景説明 + テクニックの説明」ってな感じが良いのかな
      • って、良く考えたらだいたいの人は一応そんな感じだったなw
      • 俺のLTはテクニックの説明がなかったので反省
    • 。。。と思った。id:syou6162の考えとかなりズレていることが分かりますねハイ。id:flashingwindさんの意見と若干似てるのかな。

Written by tkf

March 10, 2009 at 9:34 pm

Posted in 研究日誌

Tagged with , ,

Tsukuba.R#4で力学系の分岐現象についてLTしてきた

leave a comment »

Tsukuba.R#4でLTしてきた.

感想とか色々書きたいところだけど,疲れているので今日は資料だけUPするということで.

Rで分かる力学系 ~分岐の様子を可視化してみる~ (PDF)

Written by tkf

March 2, 2009 at 1:25 am

Posted in 研究日誌

Tagged with , ,