2015-06-26

球面三角形〜書き直し

づらづらと書いたけれど、今いちなので、その記事は下の方に回して、以下の作戦でやってみる
generating spherical eulerian triangulationというタイトルの短い論文によれば
シンプルな無向グラフの場合に限るらしいのだが
すべての球面に埋め込まれたオイラー三角化グラフは、正八面体から始めて、２つの処理を繰り返すことで作成できると言う。オイラーグラフはすべてのノードの次数が偶数であり、球面三角化の場合は面数が偶数であることにも注意する
２つの処理をｐ処理とｑ処理と呼ぶことにする
p処理
- 今、あるノードの次数が３以上であるとする(すべてのノードは偶数次数であるから実際には４以上)
- このとき、平面グラフとしてみたときに隣り合う３エッジについて次のような処理をする
- ３エッジの途中にノードを新たに作る。そのノードと残りの２エッジが連結しているノードとの間にエッジを引く
- 結果として、直線状に並んでいた３ノードが十字で仕切られたひし形に変わる
q処理
- 今、あるノードの次数が５以上であるとする
- このとき、平面グラフとしてみたときに隣り合う５得次について次のような処理をする
- 時計回りに、第３、第４番目のエッジを除去する
- その除去によって空いたスペースに第２、第５エッジが隣接するノードを結ぶエッジを引く
- その新たに引いたエッジに関して着目ノードと対側に新たにノードをおき、この新ノードに、第２、３，４，５のエッジの端点との間にエッジを引く
- 結果として、直線状に並んでいた３ノードが縦に割られた形でのひし形に変わる
Rで書いてみよう

# Spherical Eulerian Triangulation from Octahedron

library(igraph)
library(rgl)
el <- matrix(c(1,2,1,3,1,4,1,5,2,3,3,4,4,5,5,2,2,6,3,6,4,6,5,6),byrow=TRUE,ncol=2)
x <- matrix(c(0,0,1,1,0,0,0,1,0,-1,0,0,0,-1,0,0,0,-1),byrow=TRUE,ncol=3)

x <- x/sqrt(apply(x^2,1,sum))
eofv <- list()
eofv[[1]] <- c(2,3,4,5)
eofv[[2]] <- c(1,5,6,3)
eofv[[3]] <- c(1,2,6,4)
eofv[[4]] <- c(1,3,6,5)
eofv[[5]] <- c(1,4,6,2)
eofv[[6]] <- c(5,4,3,2)



my.process.choice <- function(eofv){
	current.deg <- sapply(eofv,length)
	p.cand <- which(current.deg>=4)
	q.cand <- which(current.deg>=6)
	
	s <- sample(1:length(c(p.cand,q.cand)),1)
	type <- 0
	v <- 0
	if(s <= length(p.cand)){
		type <- 1
		v <- p.cand[s]
	}else{
		type <- 2
		v <- q.cand[s - length(p.cand)]
	}
	return(list(type=type,v=v))
}

my.replace.nodes <- function(vec,v,u){
	# vecにあるvの位置をuに置き換える
	print("%%%")
	print(vec)
	print(v)
	loc <- which(vec==v)
	print("%%len")
	print(loc)
	pre <- c()
	post <- c()
	if(loc==1){
		post <- vec[-loc]
	}else if(loc==length(vec)){
		pre <- vec[-loc]
	}else{
		pre <- vec[1:(loc-1)]
		post <- vec[(loc+1):length(vec)]
	}
	c(pre,u,post)
}

rule.p <- function(eofv,v,x){
	st.v <- sample(1:length(eofv[[v]]),1)
	nb <- c(eofv[[v]],eofv[[v]])[st.v:(st.v+2)]
	print(nb)
	newv <- 1:2 + length(eofv)
	
	eofv[[nb[1]]] <- my.replace.nodes(eofv[[nb[1]]],v,c(newv[2],newv[1],v))
	eofv[[nb[2]]] <- my.replace.nodes(eofv[[nb[2]]],v,c(newv[2]))
	eofv[[nb[3]]] <- my.replace.nodes(eofv[[nb[3]]],v,c(v,newv[1],newv[2]))
	eofv[[v]] <- my.replace.nodes(eofv[[v]],nb[2],c(newv[1]))
	
	eofv[[newv[1]]] <- c(v,nb[1],newv[2],nb[3])
	eofv[[newv[2]]] <- c(newv[1],nb)
	
	newx <- rbind(x,matrix(0,2,3))
	id1 <- which(eofv[[v]]==nb[1])
	id2 <- which(eofv[[v]]==nb[2])
	id3 <- which(eofv[[v]]==nb[3])
	rest <- eofv[[v]][-(c(id1,id2,id3))]
	#newx[v,] <- (apply(matrix(x[c(id1,id2,id3),],ncol=3),2,mean) + x[v,])/2
	newx[newv[1],] <- (2*newx[v,] + newx[nb[2],])/3

	newx[newv[2],] <- (newx[v,] + 2 * newx[nb[2],])/3
	newx <- newx/sqrt(apply(newx^2,1,sum))
	return(list(eofv=eofv,x=newx))
}

rule.q <- function(eofv,v,x){
	loop <- TRUE
	while(loop){
		st.v <- sample(1:length(eofv[[v]]),1)
		nb <- c(eofv[[v]],eofv[[v]])[st.v:(st.v+4)]
		if(length(which(eofv[[nb[2]]]==nb[5])) == 0){
			loop <- FALSE
		}
		print(nb)
	}

	newv <- 1 + length(eofv)
	
	eofv[[nb[2]]] <- my.replace.nodes(eofv[[nb[2]]],v,c(newv,nb[5],v))
	eofv[[nb[3]]] <- my.replace.nodes(eofv[[nb[3]]],v,c(newv))
	eofv[[nb[4]]] <- my.replace.nodes(eofv[[nb[4]]],v,c(newv))
	eofv[[nb[5]]] <- my.replace.nodes(eofv[[nb[5]]],v,c(v,nb[2],newv))
	
	take1 <- which(eofv[[v]]==nb[3])
	take2 <- which(eofv[[v]]==nb[4])
	eofv[[v]] <- eofv[[v]][-c(take1,take2)]
	
	eofv[[newv]] <- nb[2:5]
	
	newx <- rbind(x,matrix(0,1,3))
	#rest <- eofv[[v]][-((st.v+1):(st.v+4))]
	id1 <- which(eofv[[v]]==nb[2])
	id2 <- which(eofv[[v]]==nb[3])
	id3 <- which(eofv[[v]]==nb[4])
	id4 <- which(eofv[[v]]==nb[5])
	
	rest <- eofv[[v]][-(c(id2,id3))]

	#newx[v,] <- (apply(matrix(x[rest,],ncol=3),2,mean) + x[v,])/2
	#newx[newv[1],] <- (apply(x[nb[3:4],],2,mean) + x[v,])/2
	newx[newv[1],] <- (newx[nb[3],]+newx[nb[3],])/2
	newx <- newx/sqrt(apply(newx^2,1,sum))

	return(list(eofv=eofv,x=newx))
}
el <- matrix(c(1,2,1,3,1,4,1,5,2,3,3,4,4,5,5,2,2,6,3,6,4,6,5,6),byrow=TRUE,ncol=2)
x <- matrix(c(0,0,1,1,0,0,0,1,0,-1,0,0,0,-1,0,0,0,-1),byrow=TRUE,ncol=3)

x <- x/sqrt(apply(x^2,1,sum))
eofv <- list()
eofv[[1]] <- c(2,3,4,5)
eofv[[2]] <- c(1,5,6,3)
eofv[[3]] <- c(1,2,6,4)
eofv[[4]] <- c(1,3,6,5)
eofv[[5]] <- c(1,4,6,2)
eofv[[6]] <- c(5,4,3,2)

n.step <- 1000
for(i in 1:n.step){
	proc <- my.process.choice(eofv)
	print(lapply(eofv,sort))
	print("---")
	print(proc$v)
	print(proc$type)
	if(proc$type==1){
		# rule p
		out <- rule.p(eofv,proc$v,x)
		eofv <- out[[1]]
		x <- out[[2]]
	}else{
		# rule q
		out <- rule.q(eofv,proc$v,x)
		eofv <- out[[1]]
		x <- out[[2]]
	}
}

el <- matrix(0,0,2)
for(i in 1:length(eofv)){
	el <- rbind(el,cbind(rep(i,length(eofv[[i]])),eofv[[i]]))
}
el <- unique(el)


plot3d(x,axes=FALSE)

segments3d(x[c(t(el)),])


g <- graph.edgelist(el,directed=FALSE)
degree(g)