- 最後のオイラー三角化をしながら航路を引き切る部分、力尽き果ててしまったけれど(というか、離散・乱点的アプローチより、リーマンの写像定理とメビウス変換的な共形変換ベースの航路引きの方がよさそう…と思ってしまったのでやめたのですが)
- 以下、Rmdファイル
---
title: "木とオイラー三角化のためのユーティリティ"
author: "ryamada"
date: "Sunday, May 31, 2015"
output: html_document
---
平面グラフでは、ノードにエッジが接続している順番を問題にする。それを扱うための関数をいくつか用意する。
```{r}
library(igraph)
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)
}
my.draw.graph <- function(x,el,pch=20,v.cex=1,v.col=1,seg.col=2,add=FALSE){
if(add){
segments(x[el[,1],1],x[el[,1],2],x[el[,2],1],x[el[,2],2],col=seg.col)
points(x,pch=pch,col=v.col,cex=v.cex)
}else{
plot(x,pch=pch,col=v.col,cex=v.cex)
segments(x[el[,1],1],x[el[,1],2],x[el[,2],1],x[el[,2],2],col=seg.col)
}
}
my.circle2sphere <- function(x,r=1){
len <- sqrt(apply(x^2,1,sum))
len. <- len/r*pi
xy <- sin(len.)
z <- cos(len.)
cbind(xy*x[,1]/len,xy*x[,2]/len,z)
}
my.dist.v <- function(x1,x2,type=1,r=1){
if(type==1){
return(sqrt(apply((x1-x2)^2,1,sum)))
}else if(type==2){
x1. <- my.circle2sphere(x1,r=r)
x2. <- my.circle2sphere(x2,r=r)
x1x2. <- apply(x1.*x2.,1,sum)
x1x2.[which(x1x2.>1)]<- 1
x1x2.[which(x1x2.< -1)] <- -1
return(acos(x1x2.))
}
}
```
木グラフを作る。
木グラフを作るのは簡単だが、接続具合に意識しながら、以下のように作ってみる。
```{r}
# Planar treeを作る
library(igraph)
n <- 50
el <- matrix(c(1,2),1,2)
for(i in 3:n){
el <- rbind(el,c(sample(el,1),i))
}
g <- graph.edgelist(el)
x <- layout.auto(g)
my.draw.graph(x,el)
```
頂点座標の情報を利用して、各頂点に接続しているエッジを反時計回りに並べたベクトル情報を作成する。
木が平面に埋め込まれている状態をノードIDのリスト、エッジのリスト、ノードの２次元座標として与え、
ノードに接続するエッジの順番を反時計回りに定め、その順序をエッジIDのベクトルとして返す。
```{r}
# v : node id
# el: edge list
# x : 2d coordinates
my.planar.graph <- function(v,el,x){
g <- graph.edgelist(el,directed=FALSE)
n <- length(v)
tmp.el <- t(apply(el,1,sort))
# planar graphにおいて、エッジが接続している順番を時計回りに定める
edge.order <- list()
for(i in 1:n){
tmp <- neighbors(g,i,mode=3)
tmp2 <- t(matrix(x[tmp,],ncol=2))-x[i,]
angle <- Arg(tmp2[1,] + 1i*tmp2[2,])
tmp3 <- tmp[order(angle)]
edge.order[[i]] <- tmp[order(angle)]
}
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]
}
return(list(v=v,el=el,vl=vl,x=x))
}
```
```{r}
v <- 1:n
pl.out <- my.planar.graph(v,el,x)
```
九グラフを巡回する。そのとき、あるノードに入ってきたら、進入エッジの隣のエッジ(時計回りに隣)から出ていくことにする。木グラフの場合、かならず一筆書きとなり、袋状のサイクルとなる。
このサイクルグラフを作成する。
ノードの数が２倍-1に増え、エッジの数は２倍になる。
ノードの座標は、元の木グラフの座標そのものの場合と、袋構造を反映させるために、エッジを少しずつ反時計回りに回転させて与えた直した座標の二通りを返す。
エッジの回転の強さをパラメタtheta.kで調整する。袋状であることを視認するためには少し大きな数値がよいが、そうでなければ、ごく小さな値がよい。
```{r}
my.tree2cycle <- function(pl,theta.k=1/1000){
v <- pl$v
el <- pl$el
vl <- pl$vl
x <- pl$x
e.path <- 1
e.dir.path <- 1
lb <- rep(1,length(v))
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])
}
}
x.series <- pt[2*(1:((length(pt)-1)/2))-1]
V <- 1:length(x.series) # 木を一周するノードのリスト
X <- x[x.series,] # その座標
EL <- cbind(V,c(2:length(V),1))
G <- graph.edgelist(EL)
# 座標を少しずらす
pt.mat <- matrix(pt,byrow=TRUE,ncol=2)
X.vec <- apply(X,2,diff)
x.vec.ori <- X.vec
THETA <- 2*pi/(length(x.series))*theta.k
for(i in 1:length(X.vec[,1])){
theta <- THETA*i
Rot <- matrix(c(cos(theta),-sin(theta),sin(theta),cos(theta)),byrow=TRUE,2,2)
X.vec[i,] <- Rot %*% X.vec[i,]
}
X.vec.sum <- apply(X.vec,2,cumsum)
X.vec.sum <- rbind(c(0,0),X.vec.sum)
x.vec.ori.sum <- rbind(c(0,0),apply(x.vec.ori,2,cumsum))
return(list(v=V,el=EL,x=x.vec.ori.sum,x2=X.vec.sum,x.series=x.series))
}
```
```{r}
cycle.tree.out <- my.tree2cycle(pl.out)
```
頂点座標をエッジリストを与えて、グラフを描く。
頂点座標をずらさなければ、普通の木グラフである。
```{r}
my.draw.graph(cycle.tree.out$x,cycle.tree.out$el)
```
頂点をごくわずかにずらして、袋状にしてあるが、座標のずれがごく小さいので、木に見えるのが次の描図。
```{r}
my.draw.graph(cycle.tree.out$x2,cycle.tree.out$el)
```
あえて、エッジ回転を大きめにして、袋構造を強調して描けば以下のようになる。
```{r}
cycle.tree.out2 <- my.tree2cycle(pl.out,1/10)
my.draw.graph(cycle.tree.out2$x2,cycle.tree.out2$el)
```
グラフ上の頂点からグラフに触れないようにぐるりと回って別の頂点へと曲線を引くことを考える。
木グラフからオイラー三角化を作るにはそのような曲線が引きたい。そのためにそもそもこの文書は書いている。
しかしながら、サイクルグラフがあって、その閉じた領域を陸に見立てれば、複雑な形の島の周囲を回る航路決定という意味あいに一般化できる。
以下の設定では、外周を正方形で区画し、『閉じた平面図形内の航路』という問題として設定する。
オイラー三角化復元問題の文脈では、外周の正方形は一点でに相当し、正方形が球面を開いたものとして考えていることになる。
（１）正方形内に乱点を発生させる
（２）乱点のドロネー三角化を作る
（３）ドロネー三角化の辺のうち、グラフをよぎる辺は不適格辺として削除する。グラフ非交叉三角グラフと呼ぶことにする(グラフで閉じられた内部にもサブグラフがありうる。グラフが木を袋化している場合は、そのようなものはないが、一般に、正方形内のグラフを領域の分割壁と見立てればそのようになる)
（４）グラフ非交叉三角グラフの最小全域木を作成する。これが、閉領域の『道路・航路網』とみなす。非交叉三角グラフ全体を『道路・航路網』とみなすと、陸地のごく近傍を通過する『最短沿岸路』のようなものが存在することになるが、今は、なるべく、海の真ん中を使ってほしいからである(オイラー三角化復元にあたっては、何本も航路をぶつからないように引くため)。
（５）グラフ上の２点間を結ぶ航路は、この『道路・航路網』にグラフノードを連結して作成したグラフ上の最短パスとする。ただし、グラフノードを『網』に連結するときには、ノード周囲の網点のうちグラフと非交叉のものを選ぶ。また、ノード周囲には、そのような点が選べる程度には密に乱点を発生させているものとする。
縦横の比率を変えずに１辺の長さが１の円に納まるようにグラフのノード座標を変換する
```{r}
my.standard <- function(x,scale=0.8){
ctr <- apply(x,2,mean)
x. <- t(t(x)-ctr)
len <- sqrt(apply(x.^2,1,sum))
x. <- x./max(len) * scale
x.
}
```
```{r}
scl <- 0.95
cycle.tree.out$x <- my.standard(cycle.tree.out$x,scale=scl)
cycle.tree.out$x2 <- my.standard(cycle.tree.out$x2,scale=scl)
my.draw.graph(cycle.tree.out$x,cycle.tree.out$el)
```
平面境界を構成する線分との交叉判定をするユーティリティ関数を作成する。
```{r}
# 線分の交叉判定は、ある線分が作る直線に対して残りの線分の端点が両サイドに分かれ、かつ
# もう一つの線分に対して、同じく両サイドに分かれること
my.both.side <- function(X,Y){
x1 <- X[1,1];y1 <- X[1,2];x2 <- X[2,1];y2 <- X[2,2]
X1 <- Y[1,1];Y1 <- Y[1,2];X2 <- Y[2,1];Y2 <- Y[2,2]
a1 <- (Y1-y1)*(x2-x1) - (y2-y1)*(X1-x1)
a2 <- (Y2-y1)*(x2-x1) - (y2-y1)*(X2-x1)
ret <- (a1 * a2) < 0
ret
}
my.segment.cross <- function(X,Y){
a <- my.both.side(X,Y)
b <- my.both.side(Y,X)
a & b
}
my.noncross <- function(x1,el1,x2,el2){
ret <- rep(1,length(el1[,1]))
for(i in 1:length(ret)){
X <- x1[el1[i,],]
for(j in 1:length(el2[,1])){
Y <- x2[el2[j,],]
if(my.segment.cross(X,Y)){
ret[i] <- 0
break
}
}
}
el1[which(ret==1),]
}
```
単位正方形内に乱点を発生し、そのドロネー分割をし、領域境界を定める指定グラフのエッジとの交叉を判定し、非交叉ドロネー分割を作成する。
ドロネー分割の非交叉エッジを選択する。
```{r}
library(geometry)
my.noncross.delaunay <- function(x.g,el,N=1000,Nperi=100,peri.sd = 0.2,xrange=c(0,1),yrange=c(0,1),circle=TRUE,r=1,peri.id = 1:length(x.g[,1])){
if(circle){
R1 <- matrix(runif(floor(N*4/3)*2,-r,r),ncol=2)
tmp.len <- sqrt(apply(R1^2,1,sum))
R1 <- R1[which(tmp.len <= r),]
}else{
R1 <- cbind(runif(N,xrange[1],xrange[2]),runif(N,yrange[1],yrange[2]))
}
R2 <- matrix(0,0,2)
#for(i in 1:length(x.g[,1])){
for(i in peri.id){
tmp <- cbind(rnorm(Nperi,x.g[i,1],peri.sd),rnorm(Nperi,x.g[i,2],0.1))
if(circle){
tmpid <- which(sqrt(apply(tmp^2,1,sum)) <= r)
}else{
tmpid <- which(tmp[,1]>0 & tmp[,1] < 1 & tmp[,2] > 0 & tmp[,2] < 1)
}
R2 <- rbind(R2,tmp[tmpid,])
}
R <- rbind(R1,R2)
dn <- delaunayn(R)
dn.el <- rbind(dn[,1:2],dn[,2:3],dn[,c(3,1)])
dn.el <- t(apply(dn.el,1,sort))
dn.el <- unique(dn.el)
my.draw.graph(R,dn.el)
dn.el.noncross <- my.noncross(R,dn.el,x.g,el)
return(list(x=R,el=dn.el.noncross,el.all=dn.el))
}
```
```{r}
Nrt <- 300
Nperi <- 20
peri.sd <- 0.3
dn.out <- my.noncross.delaunay(cycle.tree.out$x2,cycle.tree.out$el,Nrt,Nperi,peri.sd,circle=TRUE)
```
```{r}
my.draw.graph(dn.out$x,dn.out$el.all,seg.col="gray",v.cex=0.1)
my.draw.graph(dn.out$x,dn.out$el,seg.col=6,add=TRUE,v.cex=0.1)
my.draw.graph(cycle.tree.out$x2,cycle.tree.out$el,v.col=3,seg.col=4,add=TRUE)
```
```{r}
#edge.len <- sqrt(apply((dn.out$x[dn.out$el.all[,1],] - dn.out$x[dn.out$el.all[,2],])^2,1,sum))
edge.len <- my.dist.v(dn.out$x[dn.out$el.all[,1],],dn.out$x[dn.out$el.all[,2],])
g.ncdn <- graph.edgelist(dn.out$el.all,directed=FALSE)
edge.len <- sqrt(apply((dn.out$x[dn.out$el[,1],] - dn.out$x[dn.out$el[,2],])^2,1,sum))
g.ncdn <- graph.edgelist(dn.out$el,directed=FALSE)
mst <- minimum.spanning.tree(g.ncdn,weights=edge.len)
el.mst <- get.edgelist(mst)
my.draw.graph(dn.out$x,el.mst)
my.draw.graph(cycle.tree.out$x2,cycle.tree.out$el,v.col=3,seg.col=4,add=TRUE)
```
ただし、この場合、グラフの頂点は必ずグラフのエッジに触れているので、上で用いたmy.noncross()とは別の判定条件が必要である。
```{r}
my.find.first.step <- function(nid,g.x,g.el,R){
z <- g.x[nid,]
x1 <- rbind(R,z)
nr <- length(R[,1])
el1 <- cbind(rep(nr+1,nr),1:nr)
#nc <- rep(0,nr)
tobe.rm <- which(g.el[,1]==nid | g.el[,2]==nid)
el2 <- g.el[-tobe.rm,]
my.noncross(x1,el1,g.x,el2)[,2]
}
```
グラフのノード２つについて、
それぞれ、乱点との間に辺を引き、そのうち、グラフのエッジと交叉しないものだけを残し、乱点の非交叉最小全域木に加える。
その上で再度最小全域木を作り、その最小全域木上でのパスを２つのグラフノード間の『航路』とする。
```{r}
my.shipway <- function(two.v,g.x,g.el,R,mst.el,mst.weight,type=1){
nr <- length(R[,1])
new.vid <- nr + 1:2
first1 <- my.find.first.step(two.v[1],g.x,g.el,R)
first2 <- my.find.first.step(two.v[2],g.x,g.el,R)
new.el <- rbind(mst.el, cbind(rep(new.vid[1],length(first1)),first1),cbind(rep(new.vid[2],length(first2)),first2))
#len1 <- sqrt(apply((t(R[first1,])-g.x[two.v[1],])^2,2,sum))
#len2 <- sqrt(apply((t(R[first2,])-g.x[two.v[2],])^2,2,sum))
rep.x1 <- cbind(rep(g.x[two.v[1],1],length(first1)),rep(g.x[two.v[1],2],length(first1)))
rep.x2 <- cbind(rep(g.x[two.v[2],1],length(first2)),rep(g.x[two.v[2],2],length(first2)))
len1 <- my.dist.v(R[first1,],rep.x1,type=type)
len2 <- my.dist.v(R[first2,],rep.x2,type=type)
new.weight <- c(mst.weight,len1,len2)
new.g <- graph.edgelist(new.el,directed=FALSE)
new.mst <- minimum.spanning.tree(new.g,weights=new.weight)
new.mst.el <- get.edgelist(new.mst)
new.x <- rbind(R,g.x[two.v,])
#new.weight2 <- sqrt(apply((new.x[new.mst.el[,1],]-new.x[new.mst.el[,2],])^2,1,sum))
new.weight2 <- my.dist.v(new.x[new.mst.el[,1],],new.x[new.mst.el[,2],])
sh <- get.shortest.paths(new.mst,new.vid[1],new.vid[2],output="both",weights=new.weight2)
return(list(g = new.mst,el=new.mst.el,x=new.x,add.v=new.vid,shipway=sh))
}
```
```{r}
two.v <- sample(1:n,2)
mst.weight <- my.dist.v(dn.out$x[el.mst[,1],],dn.out$x[el.mst[,2],])
new.mst <- my.shipway(two.v,cycle.tree.out$x2,cycle.tree.out$el,dn.out$x,el.mst,mst.weight)
```
```{r}
my.draw.graph(new.mst$x,new.mst$el,v.cex=0.1)
points(new.mst$x[new.mst$add.v,],pch=20,col=4,cex=2)
my.draw.graph(cycle.tree.out$x2,cycle.tree.out$el,v.col=3,seg.col=4,add=TRUE)
my.draw.graph(new.mst$x,new.mst$el[new.mst$sh$epath[[1]],],v.col=1,seg.col=5,add=TRUE,v.cex=0.1)
points(new.mst$x[new.mst$sh$vpath[[1]],],pch=20,cex=1.3,col=5)
```
航路が引けるようになった。
円を球面に貼り付けることにする。
二次元平面のユークリッド距離を用いて最小全域木を作成してきたが、円を球面に貼り付け、その球面上の大円距離に取り換えて実施してみる。
```{r}
# type = 2 にて球面上距離を選ばせる
edge.len <- my.dist.v(dn.out$x[dn.out$el.all[,1],],dn.out$x[dn.out$el.all[,2],],type=2)
g.ncdn <- graph.edgelist(dn.out$el.all,directed=FALSE)
edge.len <- my.dist.v(dn.out$x[dn.out$el[,1],],dn.out$x[dn.out$el[,2],],type=2)
g.ncdn <- graph.edgelist(dn.out$el,directed=FALSE)
mst <- minimum.spanning.tree(g.ncdn,weights=edge.len)
el.mst <- get.edgelist(mst)
my.draw.graph(dn.out$x,el.mst)
my.draw.graph(cycle.tree.out$x2,cycle.tree.out$el,v.col=3,seg.col=4,add=TRUE)
mst.weight <- my.dist.v(dn.out$x[el.mst[,1],],dn.out$x[el.mst[,2],],type=2)
new.mst <- my.shipway(two.v,cycle.tree.out$x2,cycle.tree.out$el,dn.out$x,el.mst,mst.weight,type=2)
my.draw.graph(new.mst$x,new.mst$el,v.cex=0.1)
points(new.mst$x[new.mst$add.v,],pch=20,col=4,cex=2)
my.draw.graph(cycle.tree.out$x2,cycle.tree.out$el,v.col=3,seg.col=4,add=TRUE)
my.draw.graph(new.mst$x,new.mst$el[new.mst$sh$epath[[1]],],v.col=1,seg.col=5,add=TRUE,v.cex=0.1)
points(new.mst$x[new.mst$sh$vpath[[1]],],pch=20,cex=1.3,col=5)
```
```{r}
# cycle.tree.out 形式のオブジェクト(ノードID、エッジリスト、ノードの座標xとそれを袋対応したx2とを持つリスト)について
# 2ノードを指定(袋状のノードなので、『側』に意味がある)
# 航路を返す。航路は新たなcycle.tree.outを作る
my.update.euler <- function(cto,v,Nrt=500,Nperi=10,peri.sd=0.1,peri.id=v){
dn.out <- my.noncross.delaunay(cto$x2,cycle.tree.out$el,Nrt,Nperi,peri.sd,circle=TRUE,peri.id =peri.id)
edge.len <- my.dist.v(dn.out$x[dn.out$el[,1],],dn.out$x[dn.out$el[,2],],type=2)
g.ncdn <- graph.edgelist(dn.out$el,directed=FALSE)
mst <- minimum.spanning.tree(g.ncdn,weights=edge.len)
el.mst <- get.edgelist(mst)
mst.weight <- my.dist.v(dn.out$x[el.mst[,1],],dn.out$x[el.mst[,2],],type=2)
my.shipway(v,cycle.tree.out$x2,cycle.tree.out$el,dn.out$x,el.mst,mst.weight,type=2)
}
```
```{r}
Nrt <- 1000
Nperi <- 20
peri.sd <- 0.1
two.v <- sample(1:n,2)
new.mst <- my.update.euler(cycle.tree.out,two.v,Nrt=Nrt,Nperi=Nperi,peri.sd=peri.sd,peri.id=1:n)
my.draw.graph(new.mst$x,new.mst$el,v.cex=0.1)
points(new.mst$x[new.mst$add.v,],pch=20,col=4,cex=2)
my.draw.graph(cycle.tree.out$x2,cycle.tree.out$el,v.col=3,seg.col=4,add=TRUE)
my.draw.graph(new.mst$x,new.mst$el[new.mst$sh$epath[[1]],],v.col=1,seg.col=5,add=TRUE,v.cex=0.1)
points(new.mst$x[new.mst$sh$vpath[[1]],],pch=20,cex=1.3,col=5)
```
球面上にグラフがあって、そのグラフは球面を区画に分けるものであるときに、乱点発生とドロネー分割と最小全域木を使って、グラフ上の任意の２点間に適当な航路を定めることができた。
これを利用して、木グラフからオイラー三角化を復元することにする。
オイラー三角化が復元できるのは、Well-labeled 木グラフが平面に置かれているときである。
まず、Well-labeled木グラフを作る。
```{r}
# Well-labeled planar treeを作る
my.well.labeled.tree <- function(n,scale=0.95){
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)
# ルートノードを選ぶ
rt <- sample(1:n,1)
# ルートノードからのグラフ距離を基準に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))
g <- graph.edgelist(edge.list,directed=FALSE)
x <- layout.auto(g)
lb <- c(lb,1)
n <- n+1
# エッジの始点・終点のラベル値を格納した行列も作る
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])
# ノード番号とラベル値とを併せ持つ表示用labelを用意する
lb.2 <- c()
for(i in 1:n){
lb.2[i] <- paste(i,lb[i],sep="_")
}
v <- 1:n
pl.out <- my.planar.graph(v,edge.list,x)
cycle.tree.out <- my.tree2cycle(pl.out)
scl <- scale
cycle.tree.out$x <- my.standard(cycle.tree.out$x,scale=scl)
cycle.tree.out$x2 <- my.standard(cycle.tree.out$x2,scale=scl)
x.series <- cycle.tree.out$x.series
lb.series <- lb[x.series] # そのラベル
lb.2.series <- lb.2[x.series] # そのラベル2
return(list(v=1:(n+1),el=edge.list,x=x,lb=lb,lb.2=lb.2,edge.lb=edge.lb,pl.out=pl.out,cto=cycle.tree.out,x.series=x.series,lb.cycle=lb.series,lb.cycle.2=lb.2.series))
}
# ノード数
n <- 50
tree.lb <- my.well.labeled.tree(n)
```
```{r}
my.euler.triangulation <- function(tb){
max.lb <- max(tb$lb.cycle)
EL <- tb$cto$el
v.id <- tb$cto$el[,1]
G <- graph.edgelist(EL)
x.series = tb$x.series
lb <- tb$lb
lb.2 <- tb$lb.2
lb.series <- tb$lb.cycle
EL.lb <- cbind(lb.series[EL[,1]],lb.series[EL[,2]])
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(this.start.v==rooting.edge.target){
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))]
EL.x <- rbind(EL,c(this.start.v,this.target.v))
#if(this.start.v==2 & this.target.v==1){
# next
#}
#pre.start.v <- neighbors(G,this.start.v,mode=2)
post.start.v <- neighbors(G,this.start.v,mode=1)
pre.target.v <- neighbors(G,this.target.v,mode=2)
#post.target.v <- neighbors(G,this.target.v,mode=1)
new.start.v <- max(EL)+1
new.target.v <- max(EL)+2
v.id <- c(v.id,v.id[c(this.start.v,this.target.v)])
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] # そのラベル2
X <- tb$cto$x[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 <- new.EL
EL.lb <- cbind(lb.series[EL[,1]],lb.series[EL[,2]])
target.v <- which(lb.series==k1)
}
}
return(list(G=G,X=tb$cto$x,EL=EL,v.id=v.id))
}
```
```{r}
my.draw.eu.tri <- function(eu.out){
v.id <- eu.out$v.id
EL <- eu.out$EL
n.v <- max(v.id)
x.on.line <- (1:n.v)/n.v
x.on.line <- sort(runif(n.v))
t <- seq(from=0,to=1,length=100)*pi
r <- abs(x.on.line[1]-x.on.line[n.v])/2
ctr <- (x.on.line[1]+x.on.line[n.v])/2
x.outer <- r*cos(t) + ctr
y.outer <- r*sin(t)
plot(x.outer,y.outer,asp=TRUE,type="l")
for(i in 1:length(EL[,1])){
tmp <- v.id[EL[i,]]
tmp.ctr <- sum(x.on.line[tmp])/2
tmp.r <- abs(diff(x.on.line[tmp]))/2
tmp.x <- tmp.r*cos(t)+tmp.ctr
tmp.y <- tmp.r*sin(t)
points(tmp.x,tmp.y,type="l")
}
}
```
```{r}
# ノード数
n <- 100
tree.lb <- my.well.labeled.tree(n)
eu.out <- my.euler.triangulation(tree.lb)
my.draw.eu.tri(eu.out)
```
