多次元版 Self-avoiding path

  • Self-avoiding pathというのがある
  • 1度通過した点を2度以上訪れない軌跡
  • 酔歩的に行うことがある
  • 今、それのモデルとして、n次元空間上の乱点を三角化したグラフ上のSelf-avoiding pathをシミュレーションすることはできる
  • その軌道をk=1次元のSelf-avoiding 多様体であると考えることにする
  • さて
  • このようなランダムに広がるk次元多様体をその内部と辺縁とで考えることにする。k=1の場合の辺縁とは、パスの両端のノードのこと
  • このノードのどちらかから、エッジを選んでk=1次元多様体を伸ばす。そのときにSelf-avoidする
  • 多様体次元を上げても、内部と辺縁にするということは変わらない
  • 辺縁を構成するのが、k-1次元単体の集合になるだけ
  • Self-avoidingのルールはどうするか
  • 内部に取り込まれたk-1次元単体が接続するk次元単体は、もう選んではいけない。これがSelf-avoiding のルール
  • k=1の場合には、それは、軌跡の両端ではないノードに接続するエッジは選べないということ
  • まだバグがありそうな、お試しソース
library(geometry)
library(gtools)
library(igraph)
## Not run: 
# example delaunayn
#d = c(-1,1)
n.pt <- 100
pc = matrix(rnorm(3*n.pt),ncol=3)
tc = delaunayn(pc)

# example tetramesh
library(rgl)
clr = rep(1,3) %o% (1:nrow(tc)+1)
rgl.viewpoint(60,fov=20)
rgl.light(270,60)
tetramesh(tc,pc,alpha=0.7,col=clr)

my.make.elist <- function(tc){
	k <- ncol(tc)
	r <- nrow(tc)
	cmb <- combinations(k,2)
	v1 <- v2 <- c()
	for(i in 1:nrow(cmb)){
		v1 <- c(v1,tc[,cmb[i,1]])
		v2 <- c(v2,tc[,cmb[i,2]])
	}
	el <- cbind(v1,v2)
	el <- t(apply(el,1,sort))
	unique(el)
}

el <- my.make.elist(tc)
g <- graph.edgelist(el)
E <- as.matrix(get.adjacency(g))
E <- sign(E + t(E))
E[lower.tri(E)] <- 0
#plot(g)
#el <- get.edgelist(g)
S.list <- list()
n <- nrow(E)
S.list[[1]] <- diag(rep(1,nrow(E)))
cnt <- 1
loop <- TRUE
while(loop){
	if(length(S.list[[cnt]][,1])==0){
		break
	}
	tmp <- S.list[[cnt]] %*% E
	tmp2 <- which(tmp==cnt,arr.ind=TRUE)
	if(length(tmp2[,1])>0){
		cnt <- cnt
		S.list[[cnt+1]] <- matrix(0,length(tmp2[,1]),n)
		for(i in 1:length(tmp2[,1])){
			u <- tmp2[i,1]
			v <- tmp2[i,2]
			pre <- which(S.list[[cnt]][u,]==1)
			S.list[[cnt+1]][i,c(pre,v)] <- 1
		}
	}else{
		loop <- FALSE
	}
	cnt <- cnt+1
}

tmp <- sapply(S.list,nrow)




# Bを構成するk-1単体を含むk単体をリストアップし
# そのうち、『選んでも良いもの』を選び追加する
# 『選んでも良いもの』とは、選んだら『閉じてしまう』ものはダメ、とするのが、self-avoiding なルールだが、次元が上がると、それが単純ではない
# (1)選んだk単体がこれまでに選んだk単体ではないこと
# (2)選んだk単体が、これまでに選んだk-1単体を含まないこと ※ これが妥当そう
# (3)選んだk単体が、これまでに選んだk-1単体を含んでもよいが、穴を作らないこと(多様体としての種数を上げないこと)→結果としてこうなっていそうだが、集合知的にはこれは採用しにくい戦略かもしれない

# 簡単な方からやってみる
# S.list は一般次元三角化のk=0,1,... vs. k=0 単体の対応表のリスト
# pcは点の座標行列

# 重心に一番近い点を取り出そう
wc <- apply(pc,2,mean)
pc.wc.L <- apply((t(pc)-wc)^2,2,sum)
ctr.id <- which(pc.wc.L==min(pc.wc.L))[1]

# k次元多様体を広げる
# k=1なら線を伸ばす、k=2なら面を広げる
k <- 2


# Sをk単体 x (k-1)単体の帰属関係行列とする

tmp.k <- S.list[[k+1]]
tmp.k_ <- S.list[[k]]
tmptmp <- tmp.k %*% t(tmp.k_)
S <- sign(tmptmp == k)
S.ori <- sign(tmptmp == k)
# 変化を次のようにおさめる
# 広がりに含まれるk単体のリストWと
# その辺縁を構成するk-1単体のリストB
W <- B <- list()

# 重心付近に選んだ1点を含むk単体を一つ選び、初期状態とする
W[[1]] <- sample(which(S.list[[k+1]][,ctr.id]==1),1)
B[[1]] <- which(S[W[[1]],] == 1)
# Sから選ばれたk-1,k単体を除く
S[W[[1]],] <- 0
n.step <- 500

for(i in 1:n.step){
	print("----")
	current.B <- B[[i]]
	print(current.B)
	candidates <- which(matrix(S[,current.B],ncol=length(current.B))==1,arr.ind=TRUE)
	candidates.W <- candidates[,1]
	candidates.B <- current.B[candidates[,2]]
	print(candidates.W)
	print(candidates.B)
	selection <- sample(1:length(candidates.W),1)
	selected.W <- candidates.W[selection]
	selection2 <- which(candidates.W == selected.W)
	print(selection)
	print(selection2)
	
	print(selected.W)
	#BofnewW <- which(matrix(S[candidates.W[selection2],],nrow=length(selection2))==1,arr.ind=TRUE)
	BofnewW <- which(S[candidates.W[selection],]==1)
	print(BofnewW)
	selected.B <- candidates.B[selection2]
	print(selected.B)
	W[[i+1]] <- c(W[[i]],selected.W)
	B[[i+1]] <- unique(c(B[[i]],BofnewW))
	print("pre B")
	print(B[[i+1]])
	B[[i+1]] <- B[[i+1]][-(which(B[[i+1]] %in% selected.B))]
	print("post B")
	print(B[[i+1]])
	# 以前、Bだったけれど、内部化したprior-Bに連なるWは候補からはずす
	#S[selected.W,] <- 0
	toberemoved <- which(matrix(S[,selected.B]==1,ncol=length(selected.B)),arr.ind=TRUE)
	print(toberemoved)
	S[toberemoved[,1],] <- 0
}



plot3d(pc)
cmb <- combinations(k+1,2)
ss <- 1:length(W)
#ss <- sort(sample(1:length(W),10))
#ss <- 1:50
for(i in ss){
	for(j in 1:length(W[[i]])){
		tmp <- W[[i]][j]
		v <- which(S.list[[k+1]][tmp,]==1)
		segments3d(pc[v[t(cmb)],])
	}
}

ss <- length(W)
ss <- 20
for(i in ss){
	plot3d(pc)
	#tmp <- which(S.list[[k+1]][W[[i]],]==1,arr.ind=TRUE)
	#tmp2 <- c(t(tmp))
	#triangles3d(pc[tmp2,])
	v <- c()
	for(j in 1:length(W[[i]])){
		tmp <- W[[i]][j]
		tmp.v <- which(S.list[[k+1]][tmp,]==1)
		v <- c(v,tmp.v)
		#triangles3d(pc[v,])
	}
	triangles3d(pc[v,],alpha=0.8,col=1:length(v))
}