R で無限長ベクトル
http://www.ibm.com/developerworks/jp/linux/library/l-r3/index.html
の無限シーケンスをまねてやってみる。
# 初期値と数列を生成する関数を渡して、無限数列生成の関数を作る make_infinity <- function(init, func){ inf_vec <- init function(idx, len) { extend_by <- max(idx + len - 1 - length(inf_vec), 0) inf_vec <<- c(inf_vec, func(idx=length(inf_vec)+1, len=extend_by)) inf_vec[seq.int(idx, length.out=len)] } } # 数列を生成する関数は、ある位置 idx から長さ len の数列を生成 fib <- function(idx, len){ init <- c(-2,-1) # 以前の値を利用するときは、上位フレームからとってくる vec <- get("inf_vec", parent.frame(1))[idx + init] for(i in seq.int(3, length.out=len)) { vec <- c(vec, vec[i-2] + vec[i-1]) } vec[init] } inf_fib <- make_infinity(c(1,1), fib) inf_fib(3, 10) # [1] 2 3 5 8 13 21 34 55 89 144
漸化式ならともかく、ある idx の値を生成できる関数があるなら、
make_infinity_lite <- function(func) { inf_vec <- numeric(0) function(idx) { extend_idx <- idx[which(is.na(inf_vec[idx]))] inf_vec[extend_idx] <<- func(extend_idx) inf_vec[idx] } } inf_sin <- make_infinity_lite(sin) inf_sin(c(3,5,10)) # [1] 0.1411200 -0.9589243 -0.5440211
でよさそうやな。
この方が、必要な分しか計算しないから軽そう。
ついでに、漸化式で環境を使ってみたバージョン
関数が「参照渡し」っぽく動作する。
make_infinity_env <- function(init, func){ inf_vec <- init e <- environment() function(idx, len) { extend_by <- max(idx + len - 1 - length(inf_vec), 0) func(idx=length(inf_vec)+1, len=extend_by, env=e) inf_vec[seq.int(idx, length.out=len)] } } fib_env1 <- function(idx, len, env){ for(i in seq.int(idx, length.out=len)) { # i の値を env の環境で使えるように割り当て( assign("i_",i, env) evalq(inf_vec <- c(inf_vec, inf_vec[i_-2] + inf_vec[i_-1]), env) } # 必要なくなったオブジェクトの削除 rm(i_, envir=env) } fib_env2 <- function(idx, len, env){ # 環境のチェーンを付け替え c_env <- environment() parent.env(c_env) <- parent.env(env) parent.env(env) <- c_env for(i_ in seq.int(idx, length.out=len)) { evalq(inf_vec <- c(inf_vec, inf_vec[i_-2] + inf_vec[i_-1]), env) } # 付け替えた環境をもとに戻す parent.env(env) <- parent.env(c_env) } inf_fib_env1 <- make_infinity_env(c(1,1), fib_env1) inf_fib_env(5, 10) #[1] 5 8 13 21 34 55 89 144 233 377 inf_fib_env2 <- make_infinity_env(c(1,1), fib_env2) inf_fib_env(4, 10) #[1] 3 5 8 13 21 34 55 89 144 233
実際は、inf_vec の存在する環境で eval してるだけ。
fib_envのオブジェクトで eval 時にいるもの(ここでは i_)が少ないときは、fib_env1の方が短くてすむ。けど、多くなってくると fib_env2 の方が短くなる。
R をオブジェクト指向風に使ってみる
Rでは S4 メソッドでオブジェクト指向は行えるけど、手軽には使いにくい。
そこで、クロージャーと eval を使ってオブジェクト指向風なことをやってみる。
…というのは建前で、s「substitute と eval 使えばクロージャー内の関数呼べそう」と思ったから試してみただけ。
cls <- function(x, y) { show <- function() { cat("x:", x, ", y:", y, "\n") } add_x <- function(v) { x <<- x + v ; exe} add_y <- function(v) { y <<- y + v ; exe} a <- function(...) UseMethod("a") a.default <- function(v) { print("default") } a.matrix <- function(m) { print("matrix") } exe <- function(e){ eval(substitute(e)) } exe } #インスタンス生成 instance <- cls(3, 4) # show() メソッド呼び出し instance(show()) # x の値を返す instance(x) # x の値変更 instance(x <<- 10) # 環境経由の操作 environment(instance)$x environment(instance)$x <- 11 # メソッドを連続して呼び出し instance(add_x(3))(add_y(1))(show()) # 引数の型による多態性 m <- matrix(1:3) instance(a(m)) # a.matrix の呼び出し instance(a(1:3)) # a.default の呼び出し # Ruby の特異メソッド風に inc_x を追加 evalq(inc_x <- function(){x <<- x + 1; exe}, envir=environment(instance)) instance(inc_x())(show()) # サブクラス cls2 <- function(x,y,z) { # スーパークラスの定義をここで再定義 super <- environment(cls(x,y)) for(v in ls(super)) { val <- get(v, super) environment(val) <- environment() assign(v, val) } # サブクラスでの定義 show <- function() { cat("x:", x, ", y:", y, ", z: ", z,"\n") } # スーパークラスの同名のメソッド定義を利用するなら、こんな風にも書ける # show <- (function(){ show<-show; function(){ show(); cat("z: ", z,"\n")} })() add_z <- function(v) { z <<- z + v ; exe} exe } instance2 <- cls2(1,4,3) instance2(show()) instance2(x <<- 10) instance2(add_x(3))(add_y(1))(show())
カプセル化は protect なメソッドぐらいしかできそうにないので(めんどくさいから)やめとく。
サブクラス化の無駄が多い気がするな。再定義の部分どうにかならんかな?
super <- environment(cls(x,y)) e <- environment() parent.env(super) <- parent.env(e) parent.env(e) <- super
とやると、ある程度うまくいくけど、instance2(add_x(3))(show()) の show はスーパークラスの方が呼び出されてしまうし…。
まあ、実際に使うことはないだろうから、どうでもいいんだが。
ライフゲーム続き
改善点
- 速度向上
- インデックスのテーブル と 行列はベクトルであること を利用して outer を消した。そのおかげで apply の回数が 縦セル×横セル から 縦セル になった.。
- 状態の継続
- 前回は f() を実行するたびに初期状態から始まっていたが、2回目以降の f() は前回の終了状態から継続するようにした。
- 自動リセット
- おなじ状態を繰り返すようになったら、ランダムなセルを「生」にして変化がおきるようにした。
> summaryRprof(tmp) life_game2 <- function(x=40, y=40, a=2) { m <- matrix(sample(c(rep(0,9),1),x*y,TRUE), y, x) om <- matrix(0, y * x, a) #インデックスの行列作成 make_idx <- function(v) { rbind( c(tail(v, 1), head(v,-1)), v, c(tail(v,-1), head(v, 1)) ) } xidx <- make_idx(1:x) yidx <- make_idx(1:y) nextstep <- function(m) { livings <- t(sapply( 1:y, function(i){ # y×x の行列から 3×3x の行列を作って、9×x の行列にしたあと列ごとの和をとる row_livings <- colSums( matrix(m[yidx[,i], as.vector(xidx)],ncol=x) ) })) ifelse(livings==3, 1, ifelse(livings==4, m, 0)) } function(n=100) { for (i in 1:n) { # 過去 a 回と比較して、どれかと同じならランダムにセルを生にする if (any(apply(om == as.vector(m), 2, all))) { m <<- ( m+matrix( sample(c(rep(0,16),1),x*y,TRUE),y ) ) %% 2 } # 過去 a 回の状態を保存 om <<- cbind(as.vector(m), om[,seq(a)]) m <<- nextstep(m) image(m) } } } f <- life_game2() f() f()
スピードは約2倍になった。
> system.time(life_game()()) ユーザ システム 経過 6.75 3.26 10.57 > system.time(life_game2()()) ユーザ システム 経過 2.41 2.97 5.75
apply を呼ぶ回数を x×y 回から y 回 に減らしたから、もっと速くなると思ったけど、意外と速くならないな。
もう少し速くできないか、と思ってプロファイラーで調べてみると、
> summaryRprof() $by.self self.time self.pct total.time total.pct image.default 3.28 63.3 4.24 81.9 axis 0.46 8.9 0.50 9.7 colSums 0.26 5.0 0.64 12.4 plot.new 0.24 4.6 0.24 4.6
描画の image で 80% 以上時間を使っていた、これ以上の最適化は意味なさそうなのでこれで終わり。
ライフゲーム
最近 R を触ってなかったので、リハビリがてらにhttp://ja.wikipedia.org/wiki/ライフゲーム:セルオートマトンを書いてみる。
ライフゲームを選んだのは、「部分行列アクセスを使えば簡単かも」と思ったから。
life_game <- function(x=40, y=40) { m <- matrix(sample(c(rep(0,9),1),x*y,TRUE), y, x) nextstep <- function(m) { #上下左右をつなげる m1 <- rbind(m[nrow(m),], m, m[1,]) m1 <- cbind(m1[,ncol(m1)], m1, m1[,1]) outer(1:nrow(m), 1:ncol(m), function(vi,vj) { mapply( function(i,j) { #部分行列ごとに生きてるセル数を計算 living <- sum( m1[i:(i+2), j:(j+2)] ) if (living == 3) { 1 } else if (living == 4) { m[i,j] } else { 0 } }, vi, vj) }) } function(n=100) { for (i in 1:n) { m <- nextstep(m) image(m) } } } f <- life_game() f()
お手軽にできたけど、もとの行列の上下左右に行と列をくっつけてるあたりが、かっこ悪い。
しかも、outer のなかで apply 使ってるから計算量も多そう。
もう少し速度アップしてみるか。
function の置き換え
毎回 function と書くのは面倒くさいので、楽に書けるように置き換えてみた。
f <- function(...,b){ ab <- substitute(c(...)) if (length(ab) < 3 && !hasArg(b)) stop("2引数以上または'b'引数を指定して下さい") if (hasArg(b)) { a <- ab b <- substitute(b) }else{ len <- length(ab) a <- ab[-len] b <- ab[[len]] } da <- sub("c","",deparse(a)) db <- paste(deparse(b), collapse=";") eval.parent(parse(text=paste("function", da, db, sep=""))) }
deparse して形を整えて eval する、っていうことをやってる。
結構むだなことをしている感じ。
使い方
# 一番最後の引数が expr になる。 > square <- f(x, x^2) > square(11) [1] 121 # 名前つき引数で明示的にも指定できる # square <- f(x, b=x^2) # 引数なしの関数を作るときは、名前つき引数で指定 > f1 <- f(b=10) > f1() [1] 10 # 変数束縛もしてる > f2 <- f(x, f(y, x+y)) > f3 <- f2(10) > f3(5) [1] 15 > f3(3) [1] 13 # 複数行も"{}"で囲めばできる > f4 <- f(x,y,{ + l <- x:y + l + x + }) > f4(1,3) [1] 2 3 4 # sapply とかで使うと、こんな感じ > sapply(1:10, f(x,x+10)) [1] 11 12 13 14 15 16 17 18 19 20
感想
ちょっとした無名関数をつくるときに便利になった。
エラー処理をやってないから、関係ないエラーが出やすくなった。
ちょっと気になったこと
関数の仮引数に対して、呼び出し時の引数が少なくても問題ない。
> (function(x,y){})(10) NULL
でも、関数内で使われるとダメ。
> (function(x,y){y})(10) 以下にエラー (function(x, y) { : 引数 "y" がありませんし、省略時既定値もありません
関数実行時に引数の数をチェックしてないようだ。
仮引数の変数を評価したときに、オブジェクトが存在しなければエラーになるのかな?
ついでに、引数が多いとき。
> (function(x,y){})(10,10,10) 以下にエラー (function(x, y) { : 使われていない引数 (10)
代入先の仮引数がないからエラーになってる気がする。
R の継続
R に callCC という関数があって、
callCC package:base R Documentation Call With Current Continuation Description: A downward-only version of Scheme's call with current continuation.
と書かれていたので、試してみた。
http://www.stdio.h.kyoto-u.ac.jp/~hioki/gairon-enshuu/SchemeNotes/continuation.html を見ながらRで同じことをやってみた。
まず
(+ (* 1 2) (call/cc (lambda (c) (* 3 4)))) ;==> 14 (+ (* 1 2) (call/cc (lambda (c) (* 3 (c 4))))) ;==> 6
をほとんどそのまま書いて、
> `+`(`*`(1,2), callCC(function(c) `*`(3, 4))) [1] 14 > `+`(`*`(1,2), callCC(function(c) `*`(3, c(4)))) [1] 6
関数を明示的に呼び出さないで、
> (1 * 2) + callCC(function(c) 3 * 4) [1] 14 > (1 * 2) + callCC(function(c) 3 * c(4)) [1] 6
うまくいったので、次はこれ。
(+ (* 1 2) (call/cc (lambda (c) (set! cont c) (* 3 4)))) (cont 2) ;==> 2
> (1 * 2) + (callCC( function(c){assign("cont",c,envir=.GlobalEnv); 3 * 4})) [1] 14 > cont(2) 以下にエラー cont(2) : 戻るための関数がありません,トップレベルへジャンプします
・・・だめだ。
よくわからないのでヘルプを見直すと
「A downward-only version of Scheme's call with current continuation.」
最初は下向きってのがよく分からなかったけど、どうやら callCC を呼びたしたとこより上では継続を呼び出せないみたい。(わけが分からない文だ)
けど、それなら大域脱出はできるよね、ってことで
(define (mul-cont x) (call/cc (lambda (escape) (letrec ((mul-aux (lambda (y) (display y) (newline) (let ((result (cond ((null? y) 1) ((= (car y) 0) (escape 0)) (#t (* (car y) (mul-aux (cdr y))))))) (display (list y result)) (newline) result) ))) (mul-aux x))))) (mul-cont '(1 2 3 4)) ;(1 2 3 4) ;(2 3 4) ;(3 4) ;(4) ;() ;(() 1) ;((4) 4) ;((3 4) 12) ;((2 3 4) 24) ;((1 2 3 4) 24) ;24 (mul-cont '(1 0 3 4)) ;(1 0 3 4) ;(0 3 4) ;; carが0のときに継続を呼び出して大域脱出する. ;0
をやってみた。
mulcont <- function(x){ callCC( function(escape){ mulaux <- function(y){ cat(y, "\n"); result <- if(length(y)==0){1} else if(y[1]==0){escape(0)} else{ y[1] * mulaux(y[-1])} cat(c(y, result), "\n") result } mulaux(x) }) } > mulcont(c(1,2,3,4)) 1 2 3 4 2 3 4 3 4 4 1 4 4 3 4 12 2 3 4 24 1 2 3 4 24 [1] 24 > mulcont(c(1,0,3,4)) 1 0 3 4 0 3 4 [1] 0
よし、うまく動いた。
まだ継続自体を理解しきれてないから callCC の定義を見てもよく分からない。
> callCC function (fun) { value <- NULL delayedAssign("throw", return(value)) fun(function(v) { value <<- v throw }) } <environment: namespace:base>
いろいろ動かしながら勉強してみよ。