遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

サッカーボール：単純なルールでself-avoiding path

MIKU フラーレン サッカーボール self-avoiding path 正二十面体 黄金比
  • 昨日は平面に正六角形を敷き詰めてself-avoiding pathを描いた
  • 今日は、元の話題に戻って、サッカーボール・フラーレン上にしよう
  • サッカーボールは正２０面体のすべての面(正三角形)の角を落として正六角形を残したものなので、まずは正２０面体を作る 

# 正２０面体は１２個の頂点を持つ
# その頂点座標は以下のように簡単に出せる
G <- (1+sqrt(5))/2
g1 <- c(1,G,0)
g2 <- c(-1,G,0)
g3 <- c(1,-G,0)
g4 <- c(-1,-G,0)
X1 <- matrix(c(g1,g2,g3,g4),byrow=TRUE,ncol=3)
X2 <- X1
X2[,1] <- X1[,2]
X2[,2] <- X1[,3]
X2[,3] <- X1[,1]
X3 <- X2
X3[,1] <- X2[,2]
X3[,2] <- X2[,3]
X3[,3] <- X2[,1]

X <- rbind(X1,X2,X3)
# 一辺の長さを１に調整する
X <- X/2
# 辺となる点のペアは、全ペアのうち、長さが１(最小)のペア
D <- as.matrix(dist(X))
D[upper.tri(D)] <- 0
E <- which(D==1,arr.ind=TRUE)
plot3d(X)

for(i in 1:length(E[,1])){
	segments3d(X[E[i,],])
}
  • サッカーボール座標
    • 正２０面体のすべての辺を三等分する点を結んだのがサッカーボール座標
      • 頂点座標を求めた後、頂点ペアのうち距離が最小になるものがエッジ
      • エッジは２種類ある。６角形の２面にはさまれるエッジと６角形と５角形とにはさまれるそれ
      • ６角形にはさまれるエッジは、それが連結する頂点において、２本のエッジと120度で交叉するのに対し、６角形・５角形にはさまれるエッジは頂点での交叉角が120,108と等しくない。これを利用してエッジの色を塗り分ける

# 頂点座標を納める
Y <- matrix(0,60,3)
cnt <- 1
# 正２０面体のすべてのエッジについて三等分点を格納する
for(i in 1:length(E[,1])){
	tmp <- X[E[i,],]
	Y[cnt,] <- tmp[1,]/3 + tmp[2,]/3*2
	cnt <- cnt + 1
	Y[cnt,] <- tmp[2,]/3 + tmp[1,]/3*2
	cnt <- cnt + 1
}
# 頂点のペアワイズ距離を使って、サッカーボールのエッジを決める
D.Y <- dist(Y)
Y <- Y/min(D.Y)
D.Y <- as.matrix(dist(Y))
D.Y[upper.tri(D.Y)] <- 0
E.Y <- which(abs(D.Y-1)<10^(-3),arr.ind=TRUE)
# グラフのパッケージを使って、エッジの隣接関係を取り出す
library(igraph)
g <- graph.edgelist(E.Y)
ad.edge <- get.adjedgelist(g)
# エッジのベクトルを作る
edge.vector <- matrix(0,length(E.Y[,1]),3)
for(i in 1:length(E.Y[,1])){
	edge.vector[i,] <- apply(Y[E.Y[i,],],2,diff)
}
# エッジのタイプを分ける
e.type <- rep(1,length(E.Y[,1]))
for(i in 1:length(ad.edge)){
	tmp <- edge.vector[ad.edge[[i]],]
	tmp2 <- tmp%*%t(tmp)
	diag(tmp2) <- 0
	tmp2 <- abs(tmp2)
	for(j in 1:length(ad.edge[[i]])){
		tmp3 <- tmp2[j,which((1:3)!=j)]
		if(diff(range(tmp3))<10^(-2)){
			e.type[ad.edge[[i]][j]] <- 2
		}
	}
}


plot3d(Y)
for(i in 1:length(E.Y[,1])){
	segments3d(Y[E.Y[i,],],col=e.type[i])
}
  • まずはルールなき酔歩をしてみよう、行きつ戻りつもあり

adj.m <- as.matrix(get.adjacency(g))
adj.m <- adj.m + t(adj.m)
n.iter <- 400
v.hx <- rep(0,n.iter)
v.hx[1] <- sample(1:60,1)
v.hx[2] <- sample(which(adj.m[v.hx[1],]==1),1)


for(i in 3:n.iter){
	cand <- which(adj.m[v.hx[i-1],]==1)
	v.hx[i] <- sample(cand,1)
}
Y.walk <- Y[v.hx,]
plot3d(Y.walk)
for(i in 1:(length(Y.walk[,1])-1)){
	segments3d(Y[c(v.hx[i],v.hx[i+1]),])
  • さてルールつき酔歩
    • まだ、バグっている…
# ad.edge
# e.type
adj.m <- as.matrix(get.adjacency(g))
adj.m <- adj.m + t(adj.m)
n.iter <- 400
v.hx <- rep(0,n.iter)
v.hx[1] <- sample(1:60,1)
v.hx[2] <- sample(which(adj.m[v.hx[1],]==1),1)


for(i in 3:n.iter){
	cand <- which(adj.m[v.hx[i-1],]==1)
	v.hx[i] <- sample(cand,1)
}
Y.walk <- Y[v.hx,]
plot3d(Y.walk)
for(i in 1:(length(Y.walk[,1])-1)){
	segments3d(Y[c(v.hx[i],v.hx[i+1]),])
}

adj.m.edgeid <- adj.m
for(i in 1:length(E.Y[,1])){
	adj.m.edgeid[E.Y[i,1],E.Y[i,2]] <- i
	adj.m.edgeid[E.Y[i,2],E.Y[i,1]] <- i
	
}
n.iter <- 2000
v.hx <- rep(0,n.iter)
v.visit <- rep(0,60)
e.hx <- e.type.hx <- rep(0,n.iter-1)
v.hx[1] <- sample(1:60,1)
v.hx[2] <- sample(which(adj.m.edgeid[v.hx[1],]>0),1)
e.hx[1] <- adj.m.edgeid[v.hx[1],v.hx[2]]
e.type.hx[1] <- e.type[e.hx[1]]
v.visit[v.hx[2]] <- 1
cnt <- 3
for(i in 3:n.iter){
	print(i)
	print(v.hx[i-1])
	cand.v <- which(adj.m.edgeid[v.hx[cnt-1],]>0)
	cand.e <- adj.m.edgeid[v.hx[cnt-1],cand.v]
	cand.e.type <- e.type[cand.e]
	if(e.type.hx[cnt-2]==1){
		tmp <- which(cand.e.type==2)
		if(v.visit[cand.v[tmp]]==1){
			cnt <- cnt-2
			print("---")
			print(i)
			v.visit[v.hx[cnt-1]]<-0
			v.hx[(cnt-1)] <- 0
			e.hx[cnt-2] <- 0
			e.type.hx[cnt-2] <- 0
		}else{
			v.hx[cnt] <- cand.v[tmp]
			v.visit[v.hx[cnt]] <- 1
			e.hx[cnt-1] <- cand.e[tmp]
			e.type.hx[cnt-1] <- 2
		}
	}else{
		#cand.v2 <- which((adj.m.edgeid[v.hx[cnt-1],]>0) & (v.visit==0))
		cand.v2 <- which(adj.m.edgeid[v.hx[cnt-1],]>0)
		print(cand.v2)
		sh.cand.v2 <- sample(cand.v2)
		tmp <- sh.cand.v2[1]
		if(v.visit[tmp]==1){
			print("&&&")
			cnt <- cnt-1
			print(i)
		}else{
			v.hx[cnt] <- tmp
			v.visit[v.hx[cnt]] <- 1
			e.hx[cnt-1] <- adj.m.edgeid[v.hx[cnt-1],tmp]
			e.type.hx[cnt-1] <- 1
		}
		print(i)
	}
	cnt <- cnt+1

	if(v.hx[cnt] == v.hx[1]){
		break
	}
}
if(v.hx[cnt]==0){
	cnt <- cnt-1
}
v.hx.1 <- v.hx[1:cnt]
plot3d(Y[v.hx.1,])
for(i in 1:(length(v.hx.1)-1)){
	#if(v.hx.1[i]==0){
	#	break
	#}
	col <- e.type.hx[i]
	segments3d(Y[c(v.hx.1[i],v.hx.1[i+1]),],col=e.type.hx[i])
}