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>

いろいろ動かしながら勉強してみよ。