- 木のグラフのルートノードから出発して、常に木を右手み見るようにして木の外周を巡れば、必ずルートノードに戻ってくるし、すべてのノードでエッジのたどり方が時計回りになっている。軌跡総距離は、全エッジの長さの和の２倍である
- これを木の時計まわりの巡回とでも呼ぶことにする
- これを利用して、任意のWell-labeled 木(すべての隣接ノードが整数値のラベルを持ち、隣接ラベルの整数値の差は必ず１であるようなもの。そしてルートの値が最小(０とかにするのが普通))
- このようなWell-labeled 木について、次の規則でエッジを加えると球面上を覆うオイラー三角化が達成できることが知られている
- ラベル値が小さい方から、k+1,k+2をラベル値とするノードを両端とするエッジについて、そのエッジからグラフを時計回り巡回してラベル値がkのノードを探索する。最初に見つかったノードにエッジの両端から辺を挿入する。この挿入は平面グラフを保ったまま挿入できる
- これをコンピュータ実装しようと考えた
- ルールはきちんと決まっているので、次のような戦略をまず考えてみた
- 平面グラフであるので、グラフであるだけでなく、エッジがノードに接続しているその順番情報を付ける必要がある。平面グラフになることがわかっているグラフであるのならば、必要な情報はそれだけである。したがって、この情報を使いながら、ルールを守るように手順をコード化すればよい
- できなくはないが、面倒であることはすぐにわかった
- 何か、グラフであることを利用して、比較的単純に(エッジの接続順のような通常のグラフ情報ではないようようなものを使わずに、いわゆるグラフとしての制約だけで)実装することは可能だろうか
- それをするために考えたのが、木をぐるりとめぐるサイクルグラフに書き換えてグラフ処理をすることである
- 初めの木の状態においてはエッジの接続順の情報は必要である。それを用いて、ルートノードから出発しルートノードに戻るノード数が木の２倍-１個のグラフが作成できる
- こうすることによって、すべてのエッジは始点・終点を入れ替えた二つのエッジになる。また、その始点・終点ノードIDはそれぞれ異なり、計４個のノードによって構成される
- 一つのエッジが二つのエッジとして扱われることは、ラベル値がk+1,k+2のエッジのうち、ラベル値がkのノードへとエッジを引くべきなのは、始点のラベルがk+1,終点のラベルがk+2である。２方向のエッジが別のアイテムとして扱われるので、探索対象設定と探索処理は単純になる
- また、木グラフにエッジを加え、三角化を進めていくことで、より大きなラベル値を持つエッジに辺を加える処理のための、辺の引き先には制約が増えていく。平面グラフの用語で言えば、挿入するエッジは先行して引いてあるエッジと交叉してはいけない、ということである。その制約はグラフにグラフオブジェクトとして陽に含まれていないので、面倒くさい。その代り、辺を挿入するにあたり、グラフを分離し、複数のサイクルの集合に変換していくことで、到達可能なノードと不可能なノードとに分けることができる。これはオイラー三角化による球面の区画化そのものをグラフ的に表現したものであり、解釈も容易である
- さて、実装
- グラフ操作関数を作る
- エッジの接続順情報を用いて、木を時計回り巡回するときのノード順を作りたい
- 流入エッジの情報から、時計回り的に次に相当するエッジを流出エッジとして選択して返す関数my.next.edge()
- あるエッジからスタートし時計回りルールでエッジ的パスを作成する関数my.graph.excursion2()とそれを繰り返し使うmy.graph.excursion3()
my.next.edge <- function(el,vl,v,e,direction){
tmp <- which(vl[[v]]==e)
if(direction==1){
ret <- vl[[v]][tmp %% length(vl[[v]]) + 1]
}else{
tmp2 <- tmp - 1
if(tmp2 ==0){
tmp2 <- length(vl[[v]])
}
ret <- vl[[v]][tmp2]
}
ret.e.dir <- 1
if(el[ret,2]==v){
ret.e.dir <- 2
}
return(list(ret.e=ret,ret.e.dir=ret.e.dir))
}
my.graph.excursion <- function(el,vl,lb,e.path,e.dir.path,goal.lb,direction=1){
st.col <- 1
end.col <- 2
if(e.dir.path[length(e.dir.path)]==2){
st.col <- 2
end.col <- 1
}
st.e <- e.path[length(e.path)]
v.st <- el[st.e,st.col]
v.end <- el[st.e,end.col]
next.e <- my.next.edge(el,vl,v.end,e.path[length(e.path)],direction)
e.path <- c(e.path,next.e$ret.e)
e.dir.path <- c(e.dir.path,next.e$ret.e.dir)
new.v <- el[next.e$ret.e,3-next.e$ret.e.dir]
new.lb <- lb[new.v]
root.v <- el[e.path[1],e.dir.path[1]]
if(new.lb==goal.lb | new.v==root.v){
return(list(e.path=e.path,e.dir.path=e.dir.path))
}else{
return(my.graph.excursion(el,vl,lb,e.path,e.dir.path,goal.lb,direction=direction))
}
}
my.graph.excursion2 <- function(el,vl,lb,e.path,e.dir.path,goal.lb,direction=1){
st.col <- 1
end.col <- 2
if(e.dir.path[length(e.dir.path)]==2){
st.col <- 2
end.col <- 1
}
st.e <- e.path[length(e.path)]
v.st <- el[st.e,st.col]
v.end <- el[st.e,end.col]
next.e <- my.next.edge(el,vl,v.end,e.path[length(e.path)],direction)
e.path <- c(e.path,next.e$ret.e)
e.dir.path <- c(e.dir.path,next.e$ret.e.dir)
return(list(e.path=e.path,e.dir.path=e.dir.path))
}
my.graph.excursion3 <- function(el,vl,lb,e.path,e.dir.path,goal.lb,direction=1){
loop <- TRUE
cnt <- 1
while(loop){
out <- my.graph.excursion2(el,vl,lb,e.path,e.dir.path,goal.lb,direction)
e.path <- out$e.path
e.dir.path <- out$e.dir.path
cnt <- cnt + 1
if(cnt > 2*length(el[,1])){
loop <- FALSE
}
new.v <- el[e.path[length(e.path)],3-e.dir.path[length(e.dir.path)]]
new.lb <- lb[new.v]
if(new.lb==goal.lb){
loop <- FALSE
}
}
return(out)
}
-
- ついで、エッジパスからノードパスを作る処理と、ノードパスに登場する同一ノードに別のノードIDを割り付ける処理
- これは関数化していないので、実例ベースで以下に示す
- まず、木グラフを作り、適当にWell-labeled なラベルをつける。そのうえで、単一の"0"ノード相当のノードと、それに接続する唯一のノードがあるようにする。その心は、オイラー三角化の木グラフ表現というのは、１をラベル最小値とするWell-labeled木であって、そこから三角化を再構成するときにはラベル値１のノードのうち一つを取り立てて、そこにラベル値０のノードを接続することからスタートするのだが、その０ラベルノードを含む木を作ろうとしている。ただし、ここでは、０ラベルが１、１ラベルが２・・・としてあることに注意
- その上で、エッジパスからノードパスに変えて、巡回グラフを作成する
library(igraph)
n <- 50
edge.list <- matrix(c(1,2),1,2)
for(i in 3:n){
edge.list <- rbind(edge.list,c(sample(edge.list,1),i))
}
g <- graph.edgelist(edge.list)
plot(g)
ad.mat <- as.matrix(get.adjacency(g))
ad.mat2 <- ad.mat + (-1)*t(ad.mat)
rt <- sample(1:n,1)
lb <- shortest.paths(g,rt)+1
tmp.val <- sample(min(lb):max(lb),1)
lb <- abs(lb-tmp.val)+2
tmp.1 <- which(lb==2)[1]
edge.list <- rbind(edge.list,c(tmp.1,n+1))
lb <- c(lb,1)
n <- n+1
g <- graph.edgelist(edge.list)
plot(g)
ad.mat <- as.matrix(get.adjacency(g))
ad.mat2 <- ad.mat + (-1)*t(ad.mat)
edge.lb <- cbind(lb[edge.list[,1]],lb[edge.list[,2]])
edge.inv <- edge.lb[,1] > edge.lb[,2]
edge.list[edge.inv,] <- cbind(edge.list[edge.inv,2],edge.list[edge.inv,1])
edge.lb[edge.inv,] <- cbind(edge.lb[edge.inv,2],edge.lb[edge.inv,1])
lb.2 <- c()
for(i in 1:n){
lb.2[i] <- paste(i,lb[i],sep="_")
}
plot(g,vertex.label=lb)
pl.layout <- layout.auto(g)
plot(g, layout=pl.layout , vertex.label=lb)
edge.order <- list()
for(i in 1:n){
tmp <- neighbors(g,i,mode=3)
tmp2 <- t(matrix(pl.layout[tmp,],ncol=2))-pl.layout[i,]
angle <- Arg(tmp2[1,] + 1i*tmp2[2,])
tmp3 <- tmp[order(angle)]
edge.order[[i]] <- tmp[order(angle)]
}
plot(g, layout=pl.layout,vertex.label=lb.2)
el <- edge.list
tmp.el <- t(apply(el,1,sort))
lb <- lb
vl <- list()
for(i in 1:length(edge.order)){
vl[[i]] <- c(0)
for(j in 1:length(edge.order[[i]])){
if(i<= edge.order[[i]][j]){
vl[[i]] <- c(vl[[i]],which(tmp.el[,1]==i & tmp.el[,2]==edge.order[[i]][j]))
}else{
vl[[i]] <- c(vl[[i]],which(tmp.el[,2]==i & tmp.el[,1]==edge.order[[i]][j]))
}
}
vl[[i]] <- vl[[i]][-1]
}
e.path <- sample(1:length(el[,1]),1)
e.dir.path <- 1
goal.lb <- sample((min(lb):max(lb))[-lb[e.path[1]]],1)
goal.lb <- 0
direction <- 2
path.out <- my.graph.excursion3(el,vl,lb,e.path,e.dir.path,goal.lb,direction)
pt <- c()
for(i in 1:length(path.out[[1]])){
tmp <- el[path.out[[1]][i],]
if(path.out[[2]][i]==1){
pt <- c(pt,tmp)
}else{
pt <- c(pt,tmp[2:1])
}
}
pt
x.series <- pt[2*(1:((length(pt)-1)/2))-1]
v <- 1:length(x.series)
lb.series <- lb[x.series]
lb.2.series <- lb.2[x.series]
X <- pl.layout[x.series,]
X.ori <- X
EL <- cbind(v,c(2:length(v),1))
G <- graph.edgelist(EL)
EL.lb <- cbind(lb.series[EL[,1]],lb.series[EL[,2]])
plot(G,layout=X.ori+rnorm(length(X.ori),0,0.1),vertex.size=3)
plot(G, layout=X.ori+rnorm(length(X.ori),0,0.1),vertex.label=lb.2.series)
max.lb <- max(lb)
for(i in 3:max.lb){
k1 <- i-2
k2 <- i-1
k3 <- i
select.e <- which(EL.lb[,1]==k2 & EL.lb[,2]==k3)
start.v <- c(t(EL[select.e,c(2,1)]))
target.v <- which(lb.series==k1)
for(j in 1:length(start.v)){
if(i==3){
if(j==2){
next
}
}
this.start.v <- start.v[j]
tmp <- shortest.paths(G,this.start.v,target.v)
this.target.v <- target.v[which(tmp==min(tmp))]
post.start.v <- neighbors(G,this.start.v,mode=1)
pre.target.v <- neighbors(G,this.target.v,mode=2)
new.start.v <- max(EL)+1
new.target.v <- max(EL)+2
x.series <- c(x.series,x.series[c(this.start.v,this.target.v)])
lb.series <- lb[x.series]
lb.2.series <- lb.2[x.series]
X <- pl.layout[x.series,]
new.edges <- rbind(c(this.start.v,this.target.v),c(pre.target.v,new.target.v),c(new.target.v,new.start.v),c(new.start.v,post.start.v))
new.EL <- rbind(EL,new.edges)
rm.edge <- c(which(new.EL[,1]==this.start.v & new.EL[,2]==post.start.v),which(new.EL[,1]==pre.target.v & new.EL[,2]==this.target.v))
print(rm.edge)
print(dim(new.EL))
new.EL <- new.EL[-rm.edge,]
print(dim(new.EL))
G <- graph.edgelist(new.EL)
EL <- get.edgelist(G)
EL.lb <- cbind(lb.series[EL[,1]],lb.series[EL[,2]])
target.v <- which(lb.series==k1)
}
}
plot(G,X)
EL2 <- cbind(x.series[EL[,1]],x.series[EL[,2]])
G2 <- graph.edgelist(EL2)
plot(G)
plot(G2)
