Hatena::ブログ(Diary)

Great oaks from little acorns grow.

2009-10-07

R で無限長ベクトル

| 22:51 |  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 の方が短くなる。

2009-05-26

R をオブジェクト指向風に使ってみる

| 00:03 |  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 はスーパークラスの方が呼び出されてしまうし…。

まあ、実際に使うことはないだろうから、どうでもいいんだが。

2009-05-13

ライフゲーム続き

| 22:51 |  ライフゲーム続きのブックマークコメント

改善点

速度向上
インデックスのテーブル と 行列はベクトルであること を利用して 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% 以上時間を使っていた、これ以上の最適化は意味なさそうなのでこれで終わり。

2009-05-12

ライフゲーム

| 22:51 |  ライフゲームのブックマークコメント

最近 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 使ってるから計算量も多そう。

もう少し速度アップしてみるか。

2008-07-20

function の置き換え

| 22:51 |  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

感想

ちょっとした無名関数をつくるときに便利になった。

エラー処理をやってないから、関係ないエラーが出やすくなった。