ryamadaのコンピュータ・数学メモ RSSフィード

統計学・遺伝学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
数学用語集(TOMAC)

2017-09-13 球面オイラー三角化と Steiner triple trade

[][][]球面オイラー三角化とSteiner triple trade 08:39

  • [:title=こちら]に、球面オイラー三角化とSteiner triple tradeとの関係が書いてある
  • 簡単に言うと:
    • 三角形の頂点ノードを三つ組み(x,y,z)とする
    • オイラー三角化では、三角形が2色に塗り分けられるが、各色ごとに三角をグループ分けすることにする
    • 今、色1のグループに(x,y,z)があったとき、それは色2のグループに(x',y,z),(x,y',z),(x,y,z')が必ず帰属していることを意味する。ただし、xとx'は違う(y,zも同様)
  • このような三つ組みのペアが、オイラー三角化に対応していることを利用して、オイラー三角化を組み合わせの論理で書き進めているのが、リンク先のペイパー
トラックバック - http://d.hatena.ne.jp/ryamada/20170913

2017-09-02 Monge property

[][][]Monge property 09:01

  • 平面グラフとは、エッジを交叉させずに平面に描けるグラフのこと
  • 平面グラフの距離行列はMonge propertyという性質を満足する
  • 平面グラフの頂点u,vから、別の頂点x,yへの最短距離を考えるとき d(u,x)+d(v,y)とd(u,y)+d(v,x)とのどちらが小さいかを問題にしたい
  • 今、d(u,x)+d(v,y) <= d(u,y)+d(v,x)となっていたときに、u,vとx,yに、順序を定めて u<=v & x<=yとすることにする
  • 今、u,vは2ノード、x,yも2ノードだったが、これを、平面グラフのノードの部分集合と別の部分集合とにして、n個のノード、m個のノードの部分集合でかんがえたときにも、n個、m個内の順序関係が定まる
  • この性質がMonge property
  • この性質が成り立っているとき、巡回セールスマン問題はこの性質の成立を前提としないそれよりも簡単だという
  • また、Monge propertyの対称正方行列版をSupnick matrixと呼ぶ
  • planar graphの距離行列にこれがあるらしいのだが、部分集合の取り方の定義のところがまだ良くわかっていない…。任意の部分集合2つを取ってよいわけではないようだ(それを確かめたのが後述のコード)
  • 参考1
  • 参考2
  • 参考3
  • 参考4
library(igraph)

el <- matrix(c(1,2,2,3,1,4,1,7,2,7,2,8,3,8,7,9,7,8,8,9,3,6,9,4,9,5,9,6,8,6,4,5,5,6),byrow=TRUE,ncol=2)
g <- graph.edgelist(el,directed=FALSE)
plot(g)

n.iter <- 100
ret <- rep(0,n.iter)
ret2 <- ret
for(ii in 1:n.iter){
	L <- c(1,2,2,1,3,1,5,1,3,5,3,4,2,4,2,8,2)
	#L <- runif(17)
	d <- distances(g,weights=L)
	d. <- d[c(1,2,3),c(4,5,6)]

	tmp <- c()
	tmp2 <- tmp
	for(i in 1:2){
		for(j in (i+1):3){
			for(k in 1:2){
				for(l in (k+1):3){
					d.. <- d.[c(i,j),c(k,l)]
					tmp <- c(tmp,d..[1,1] + d..[2,2] - (d..[1,2]+d..[2,1]))
					# 行列の要素xをexp(x)にしてその行列式をとると必ず1以下
					tmp2 <- c(tmp2,det(exp(d..)))
				}
			}
		}
	}

	ret[ii] <- (prod(tmp<=0))
	ret2[ii] <- prod(tmp2<=1)
}

plot(ret)
plot(ret2)
AB.0 <- matrix(sample(1:1000,6),ncol=2)
ord <- rbind(c(1,2,3),c(1,3,2),c(2,1,3),c(2,3,1),c(3,1,2),c(3,2,1))
n <- length(ord[,1])

for(k in 1:n){
	for(kk in 1:n){
		AB <- cbind(AB.0[ord[k,],1],AB.0[ord[kk,],2])
		abcd <- matrix(0,0,4)
		val <- c()
		for(i in 1:2){
			for(ii in (i+1):3){
				for(j in 1:2){
					for(jj in (j+1):3){
						tmp <- d.[c(AB[i,1],AB[ii,1]),c(AB[j,2],AB[jj,2])]
						abcd <- rbind(abcd,c(i,ii,j,jj))
						#abcd <- rbind(abcd,c(AB[i,1],AB[ii,1],AB[j,2],AB[jj,2]))
						val <- c(val,tmp[1,1] + tmp[2,2] - (tmp[1,2]+tmp[2,1]))
					}
				}
			}
		}
		print(prod(sign(val)==1))
		if(prod(sign(val)==1)){
			print(val)
		}
	}
}


#cbind(abcd,val)

library(gtools)
n <- 4
vs <- sample(1:1000,n)
pm <- permutations(n,n)

pm2 <- permutations(n,2)
pm2 <- t(apply(pm2,1,sort))


for(i in 1:length(pm[,1])){
	tmp <- d.[vs[pm[i,]],vs[pm[i,]]]
	ret <- c()
	for(j in 1:length(pm2[,1])){
		for(jj in 1:length(pm2[,1])){
			tmp2 <- tmp[pm2[j,],pm2[jj,]]
			tmp3 <- tmp2[1,1] + tmp2[2,2] - (tmp2[1,2]+tmp2[2,1])
			ret <- c(ret,tmp3)
		}
	}
	print(prod(sign(ret)==1))
		if(prod(sign(ret)==1)){
			print(ret)
		}
}


トラックバック - http://d.hatena.ne.jp/ryamada/20170902

2017-08-31 交点

[] 16:22

  • 同一平面上にある4点が、一般的な位置関係にあるとき、4点を2点ずつのペアに分けて、ペアを通る直線を2本引いたときの交点を出したい
  • 散らし書きコード
my.2d.intersect <- function(v1,v2,v3,v4){
	M <- matrix(c(v2[1]-v1[1],v3[1]-v4[1],v2[2]-v1[2],v3[2]-v4[2]),byrow=TRUE,2,2)
	b <- matrix(c(-v1[1]+v3[1],-v1[2]+v3[2]),ncol=1)
	a <- solve(M,b)
	return(a)
}

my.3d.intersect.tri <- function(v1,v2,v3,v4){
	a <- my.2d.intersect(v1[1:2],v2[1:2],v3[1:2],v4[1:2])
	t <- a[1]
	v1 + (v2-v1) * t
}

v1 <- rnorm(3)
v2 <- rnorm(3)
v3 <- rnorm(3)
v4 <- rnorm(3)
v1 <- v1/sum(v1)
v2 <- v2/sum(v2)
v3 <- v3/sum(v3)
v4 <- v4/sum(v4)

v5 <- my.3d.intersect.tri(v1,v2,v3,v4)
plot3d(rbind(v1,v2,v3,v4,v5))
spheres3d(rbind(v1,v2,v3,v4),radius=0.1)
spheres3d(v5,radius=0.1,col=2)
トラックバック - http://d.hatena.ne.jp/ryamada/20170831

2017-08-30 三角形版「整数分割」

[][]三角形版「整数分割」 18:07

  • 整数をいくつかの整数に分けることを整数分割と言う
  • 今、辺の長さが整数であるような三角形を「整数三角形」とでも呼ぶことにする
  • ただし、ここで言う「辺の長さ」は本当の辺の長さではなく、「この辺は整数●単位分に相当する」と考えたいときの●が整数である、という意味とする
  • これをあるルールで整数分割することを再帰的に繰り返し、より小さな整数三角形でのタイリングにすることにする
  • ひとまず走り書きコード
# 三角形を(x,y,z)分割する
library(partitions)
library(combinat)
library(rgl)
# 同一平面上にある4点が作る直線の交点座標を出す
my.2d.intersect <- function(v1,v2,v3,v4){
	M <- matrix(c(v2[1]-v1[1],v3[1]-v4[1],v2[2]-v1[2],v3[2]-v4[2]),byrow=TRUE,2,2)
	b <- matrix(c(-v1[1]+v3[1],-v1[2]+v3[2]),ncol=1)
	a <- solve(M,b)
	return(a)
}

my.3d.intersect.tri <- function(v1,v2,v3,v4){
	a <- my.2d.intersect(v1[1:2],v2[1:2],v3[1:2],v4[1:2])
	t <- a[1]
	v1 + (v2-v1) * t
}

my.divide.n <- function(n){
	if(n==1){
		return(list(composition=diag(rep(1,3)),tris=matrix(c(1,2,3),nrow=3)))
	}
	comp <- compositions(n,3)
	nv <- length(comp[1,])
	trios <- matrix(combn(1:nv,3),nrow=3)
	e1 <- comp[,trios[1,]] - comp[,trios[2,]]
	e2 <- comp[,trios[2,]] - comp[,trios[3,]]
	e3 <- comp[,trios[3,]] - comp[,trios[1,]]
	tris <- which((apply(abs(e1),2,sum) == 2) & (apply(abs(e2),2,sum) == 2) & (apply(abs(e3),2,sum) == 2))
	return(list(composition=comp,tris=trios[,tris])) # n^2個の三角形ができる
}

my.divide.n(2)


my.integer.division <- function(n,m){
	p <- min(n,m)
	q <- max(n,m)
	tmp1 <- q %% p
	tmp2 <- q %/% p
	ret <- rep(tmp2,p)
	tmp3 <- floor(((1:p) * tmp1)/p)
	tmp4 <- c(0,diff(tmp3))
	tmp5 <- which(tmp4==1)
	ret[tmp5] <- ret[tmp5] + 1
	return (ret)
}

my.integer.division(3,22)

my.div.series <- function(d){
	ret <- list()
	ret[[1]] <- d
	for(i in 2:length(d)){
		n <- length(ret[[i-1]])
		tmp1 <- ret[[i-1]][1:(n-1)]
		tmp2 <- ret[[i-1]][2:n]
		ret[[i]] <- apply(cbind(tmp1,tmp2),1,min)
	}
	return(ret)
}
id <- my.integer.division(3,22)
id
my.div.series(id)

my.divide.tri <- function(tri,d.n=NULL){
	if(sum(abs(tri[[1]]-c(1,1,1)))==0){
		ret <- list()
		#ret[[1]] <- list()
		#ret[[1]][[1]] <- tri[[1]]
		#ret[[1]][[2]] <- tri[[2]]
		ret[[1]] <- list(tri[[1]],tri[[2]])
		return(ret)
	}
	if(min(tri[[1]])==1){
		snd <- sort(tri[[1]])[2]
		if(snd == 1){
			ret <- my.divide.tri.11(tri)
			return(ret)
		}else{
			ret <- my.divide.tri.1(tri)
			return(ret)
		}
	}
	
	# 等分割の場合の、三角形のリストと、そのBarycentric coordinates相当整数値
	
	if(is.null(d.n)){
		
		d.n <- my.divide.n(min(tri[[1]]))
	}
	
	# 整数分割
	edge.div <- list()
	edge.frac <- list()
	div.series <- list()
	for(i in 1:3){
		edge.div[[i]] <- my.integer.division(min(tri[[1]]),tri[[1]][i])
		div.series[[i]] <- my.div.series(edge.div[[i]])
		edge.frac[[i]] <- cumsum(c(0,edge.div[[i]]))
		edge.frac[[i]] <- edge.frac[[i]]/sum(edge.div[[i]])
	}
	

	
	# edge上の点座標の計算
	nv <- max(d.n$tris)
	divnum <- max(d.n$composition)
	v <- c(1:3,1:3,1)
	edge.x <- list()
	for(i in 1:3){
		edge.x[[i]] <- matrix(0,3,divnum+1)
		for(j in 1:(divnum+1)){
			edge.x[[i]][,j] <- (1-edge.frac[[i]][j]) * tri[[2]][,v[i+1]] + edge.frac[[i]][j] * tri[[2]][,v[i+2]]
			print(edge.x[[i]][,j])
		}
	}
	# 分割後の頂点座標の計算

	x <- matrix(0,3,nv)
	v.x <- list()
	
	vv <- matrix(0,0,3)
	for(i in 1:nv){
		tmp <- d.n$composition[,i]
		num.zero <- length(which(tmp==0))
		if(num.zero==2){
			nonzeroid <- which(d.n$composition[,i]!=0)
			x[,i] <- tri[[2]][,nonzeroid]
						print("$$$")
			print(x[,i])
		}else if(num.zero==1){
			zeroid <- which(d.n$composition[,i]==0)
			val <- tmp[v[zeroid+2]]
			x[,i] <- edge.x[[zeroid]][,val+1]
			print("&&&")
			print(x[,i])
		}else{
			v.x[[i]] <- list()
			three.lines <- list()
			for(j in 1:3){
				one <- edge.x[[v[j+1]]][,tmp[j]+1]
				another <- edge.x[[v[j+2]]][,divnum-tmp[j]+1]
				#print(one)
				#print(another)
				three.lines[[j]] <- cbind(one,another)
			}
			tmp <- cbind(three.lines[[1]],three.lines[[2]],three.lines[[3]])
			cross1 <- my.3d.intersect.tri(three.lines[[1]][,1],three.lines[[1]][,2],three.lines[[2]][,1],three.lines[[2]][,2])
			cross2 <- my.3d.intersect.tri(three.lines[[2]][,1],three.lines[[2]][,2],three.lines[[3]][,1],three.lines[[3]][,2])
			cross3 <- my.3d.intersect.tri(three.lines[[3]][,1],three.lines[[3]][,2],three.lines[[1]][,1],three.lines[[1]][,2])

			x[,i] <- apply(cbind(cross1,cross2,cross3),1,mean)

		}


		
		
	}
	x
	plot3d(t(x))
	spheres3d(t(x),radius=0.1)

	ret <- list()
	for(i in 1:length(d.n$tris[1,])){
		# 3頂点のID
		vid <- d.n$tris[,i]
		# 3頂点の座標
		vx <- x[,vid]
		
		# 3辺の長さ
		# 第1辺(第1頂点の対辺)
		tmp <- d.n$composition[,vid]
		tmptmp1 <- tmp[,1]-tmp[,2]
		tmptmp2 <- tmp[,2]-tmp[,3]
		tmptmp3 <- tmp[,3]-tmp[,1]
		
		dir1 <- which(tmptmp1==0)
		level1 <- tmptmp1[dir1]
		mm1 <- tmptmp1[v[dir1+2]]
		n1 <- div.series[[dir1]][[level1+1]][mm1+1]
		dir2 <- which(tmptmp2==0)
		level2 <- tmptmp2[dir2]
		mm2 <- tmptmp2[v[dir2+2]]
		n2 <- div.series[[dir2]][[level2+1]][mm2+1]
		dir3 <- which(tmptmp3==0)
		level3 <- tmptmp3[dir3]
		mm3 <- tmptmp3[v[dir3+2]]
		n3 <- div.series[[dir3]][[level3+1]][mm3+1]

		ret[[i]] <- list(c(n1,n2,n3),vx)
	}
	return(ret)
}
my.divide.tri.11 <- function(tri){ # tri[[1]] = {1,1,n!=1}
	ord <- order(tri)
	snd <- tri[ord[2]]
	mx <- tri[ord[3]]
	
	ret <- list()
	
	n.triangle <- mx
	for(i in 1:n.triangle){
		
		ret[[i]] <- list(c(1,1,1),cbind(x1,x2,x3))
	}
	return(ret)
}

my.divide.tri.1 <- function(tri){ # tri[[1]] = {1,n!=1,m!=1}
	ord <- order(tri)
	snd <- tri[ord[2]]
	mx <- tri[ord[3]]
	
	n.d <- my.integer.division(snd,mx)
	
	
	
}

my.plot.divtri <- function(q,face=TRUE,edge=TRUE){
	mat <- matrix(0,0,3)
	for(i in 1:length(q)){
		mat <- rbind(mat,t(q[[i]][[2]]))
	}
	plot3d(mat)
	if(face){
		col <- matrix(0,0,3)
		for(i in 1:length(q)){
			col <- rbind(col,q[[i]][[1]])
		}
		col <- col
		col <- t(t(col)/apply(col,2,max))
		col <- (1-col)

		for(i in 1:length(q)){
			#triangles3d(t(q[[i]][[2]]),col=rgb(col[i,1],col[i,2],col[i,3]))
			triangles3d(t(q[[i]][[2]]),col=i)
		}
	}

	if(edge){
		ed <- matrix(0,0,3)
		for(i in 1:length(q)){
			#segments3d(rbind(q[[i]][[2]][,1],q[[i]][[2]][,2]))
			#segments3d(rbind(q[[i]][[2]][,2],q[[i]][[2]][,3]))
			#segments3d(rbind(q[[i]][[2]][,3],q[[i]][[2]][,1]))
			ed <- rbind(ed,rbind(q[[i]][[2]][,1],q[[i]][[2]][,2]))
			ed <- rbind(ed,rbind(q[[i]][[2]][,2],q[[i]][[2]][,3]))
			ed <- rbind(ed,rbind(q[[i]][[2]][,3],q[[i]][[2]][,1]))
		}
		segments3d(ed)
	}
}

my.divide.tri.recursive <- function(tri){
	loop <- TRUE
	cnt <- 1
	
	while(loop){
		loop <- FALSE
		ret <- list()
		print("cnt")
		print(cnt)
		cnt <- cnt+1
		for(i in 1:length(tri)){
			if(sum(abs(tri[[i]][[1]]-c(1,1,1)))==0){
				cnt2 <- length(ret)+1
				#ret[[cnt2]] <- list()
				#ret[[cnt2]][[1]] <- tri[[i]][[1]]
				#ret[[cnt2]][[2]] <- tri[[i]][[2]]
				#ret[[cnt2]] <- list(tri[[i]][[1]],tri[[i]][[2]])
				ret[[cnt2]] <- tri[[i]]
				print("in")
				print(abs(tri[[i]][[1]]-c(1,1,1)))

			}else{
				print(tri[[i]])
				out <- my.divide.tri(tri[[i]])
				print(out)
				for(j in 1:length(out)){
					cnt2 <- length(ret)+1
					#ret[[cnt2]] <- list()
					#ret[[cnt2]][[1]] <- out[[j]][[1]]
					#ret[[cnt2]][[2]] <- out[[j]][[2]]
					#ret[[cnt2]] <- list(out[[j]][[1]],out[[j]][[2]])
					ret[[cnt2]] <- out[[j]]
				}
				#ret <- my.divide.tri.recursive(tri=out,ret=ret)
				#print(j)
				print(length(ret))
				print("----")
				print(ret)
				print("===")
				#loop <- TRUE
			}
		}
		#print(tri)
		tri <- ret
		#print(tri)
	}

	return(tri)
}

tri.org <- list(le = c(5,7,8),x = diag(rep(3,3)))
tri <- list()
tri[[1]] <- list(tri.org[[1]],tri.org[[2]])

#tri[[1]] <- list()
#tri[[1]][[1]] <- tri.org[[1]]
#tri[[1]][[2]] <- tri.org[[2]]


tmpout <- my.divide.tri(tri[[1]])

my.plot.divtri(tmpout)
for(i in 1:length(tmpout)){
	ctr <- apply(tmpout[[i]][[2]],1,mean)
	txt = paste("",i)
	text3d(ctr+0.1,text=txt)
}



tmpout2 <- my.divide.tri.recursive(tri)

my.plot.divtri(tmpout2)
トラックバック - http://d.hatena.ne.jp/ryamada/20170830

2017-08-26 『私のためのBrownian map構成法』

[][][][][]p-angulation n faces平面グラフとそのグラフ距離『私のためのBrownian map構成法』 16:54

  • Uniqueness and universality of the Brownian map(資料1)
  • Random Geometry on the Sphere(資料2)
  • The Brownian map is the scaling limit of uniform random plane quadrangulations(資料3)
  • 平面グラフは、球面上の連結グラフであるので、平面グラフはorientation-preserving homeomorphisms of S2と考えられる
  • 平面グラフには双対グラフがある(多角形をノードにとり、多角形の隣接関係をエッジとしたグラフ)
  • 平面グラフには次のような制約を入れられる
    • p角形のみで出来ていて、n個のp角形で出来ているもの
    • この集合をA^q_nと書く。A set of all rooted planar q-angulations with n faces と英語で記載するとわかりやすいかもしれない
    • q=3のとき、nが奇数のときはA^{q=3}_nには要素がない
    • 集合A^p_nの要素M_n ¥in  A^q_n
  • M_nはグラフなので、ノードのセットがある。それをm_nと書き、そのグラフのノード間距離(すべてのエッジの長さは1)をd_{gr}と書けば、それはこのグラフを「表している」から、(m_n,d_{gr})と書くことにする。
    • グラフがG=(V,E)のように、ノード集合とエッジ集合とのペアであることと対応するが、エッジ集合の代わりにd_{gr}を用いている。
  • グラフ同士に「異同」を定量する『メトリック』がある
    • 2つの(m_n^1,d_{gr}^1),(m_n^2,d_{gr}^2)があったときに、それらが全く同一のグラフであったときに0となり、異同の程度に応じて正の値をとるような、そんなメトリックがある。Gromov-Hausdorff dissimilarityと言うのがそれである
    • このとき、(m_n,d_{gr})は、ある空間Kの点であって、この空間KにはメトリックD_{GH}が備わっている、と言い、(K,d_{GH})と書いて、そのゆおなメトリックを備えた空間であることを表記する
  • q-angulationでのfaces数nをどんどん増やしていくと、対応する球面は、点の粗密はあるものの、点で埋め尽くされるようになる。そのとき、どの点とどの点とが近くて、どの点とどの点とが遠いか、というのは、グラフ距離d_{gr}で測れるのは、相変わらずである。ただし、点の数がどんどん増えると、グラフ的距離はどんどん長くなる。したがって、nの値に応じて、1辺の長さを加減するのが適当である。
    • 今、(¥frac{9}{q(q-1)})^{1/4} n^{-1/4}という調整係数を用いることにすると
    • (m_n, (¥frac{9}{q(q-1)})^{1/4} n^{-1/4} d_{gr})というグラフ(とその遠近メトリックと)が定められるが、このnを無限大にしたときの極限が、qによらず同じであることが知られているのだという
    • この極限をブラウニアンマップと呼ぶ
    • それを(m_n,(¥frac{9}{q(q-1)})^{1/4} n^{-1/4} d_gr) ¥rightarrow_{n ¥to ¥infty}(m_¥infty,D^*)と書く
    • ここで、nが増えれば増えるほど、平均して、点間距離はn^{-1/4}に比例して長さを調整する必要がある、ということが読み取れる。単純に2次元平面なら、-1/2乗に比例させるのがよいが、グラフ距離の場合には、次元が2ではなく、もっとフラクタル的な込み入った次元なので-1/4となる(らしい)

[][][][][]酔歩の掛け合わせとしてのContinuous Random tree とブラウニアンマップ『私のためのBrownian map構成法』 16:54

  • ブラウン散歩を木とみなすことで出来るのがCRT(Continuum Random Tree)
  • そのCRTの上に酔歩を乗せて、それによって、木の各所に値をラベル付けすることで、「ラベル付けされたCRT」が出来上がる
  • この「ラベル付けされたCRT」の点集合に点間距離が定義づけできて、それがブラウニアンマップである
  • 次の2つのサブ記事で、その2点を解説する

[][][][][]ブラウン散歩はCRT『私のためのBrownian map構成法』 17:04

  • ブラウン散歩は、0から出発して0に戻る酔歩であって、正の値のみをとっているもの
  • 特に、0 ¥le s ¥le 1の時間範囲に限定したものがNormalized Brownian Excursion(標準ブラウン散歩)
  • 今、時刻s,tにおける位置をe_s,e_tのように書くことにしとする
  • d_e(s,t) = e_s + e_t - 2 min_{min(s,t) ¥le r ¥max(s,t)}(e_r)という値を、時刻sと時刻tとの「ブラウン散歩eの意味での距離」と定義する
    • s=tのときd_e(s,t=s)=0だし、
    • s,tの間の値が、すべてe_s,e_t以上である場合も、d_e(s,t)=0であるから
    • ブラウン散歩で原点に近づく向きに動いている経路は、必ず、原点から遠ざかる経路とペアになってd_e(s,t)=0なる点とペアになっている
  • 結局、原点0から出発し、原点0に戻るような「木」が得られる
  • 円周の弧をあちこちでつぶして木を作ったようなものでもある
  • このつぶした円周は、「円だったとしたら時計回りに回る」ような動き方を「木の上」で行うことができて、それは、木の外側をなぞりながら木を一周するものに相当する
  • T_e: S^1/¥sim _eと書いて、単位円S^1を、ブラウン散歩eという意味での点の同一視関係¥sim_eによって商を取ったものが、CRT T_e(eというブラウン散歩によるCRT木である、と読む
  • ブラウン散歩では、時刻パラメタ0 ¥le s ¥le 1があった
  • T_eは木であるので、その上の点a ¥in T_eのように書ける
  • ブラウン散歩eとそれに対応するCRTT_eとの間には、次のような対応付けマッピングがあるp_e(u) = a; 0 ¥le u ¥le 1, a ¥in T_e
    • 特に、木の上の一つの点a ¥in T_eにはp_e(u_1) = a, p_e(u_2) = aのように、複数の0 ¥le u_i ¥le 1が対応づく。さらに特に、T_e上で枝分かれするところでは3個(以上)の0 ¥le u_i ¥le 1が対応する
  • 木の周回弧とそれに対応する、ブラウン散歩の時間セグメントの関係は次のように定める
    • 木の上の二点a,b ¥in T_eがあるとき、aからbへと木を時計回りに周回することを考える。周回なので、木の枝のこちら側かあちら側かの区別ができるので、複数の周回の仕方がある
    • それはp(u_a1)=p(u_a2)=ap(u_b1)=p(u_b2)=bとで考えれば、u_a1-> u_b1, u_a1-> u_b2,u_a2 -> u_b1, u_a2 -> u_b2の4通りがあるようなものである。この4通りについて、0 ¥le u_*i ¥le 1なる時刻パラメタで考えたときの最短時間を、T_e上の2点の時計回りの弧と定め、¥[a,b¥]と書くことにする。
    • ¥[b,a¥]と木の上の2点を入れ替えると、ルートノード(s=0に対応する点)を通過するような周回路が対応する
    • この¥[a,b¥]は、次の記事での、ラベル付けされたCRT上の2点間の距離の定義で使用する

[][][][][]CRT上の酔歩ラベル付けとそれが表すメトリック〜ブラウニアンマップ『私のためのBrownian map構成法』 17:04

  • ブラウニアンラベルでは、T_e上に値を貼り付ける
  • ブラウン散歩eでの時刻0での点に対応するT_eの点¥rhoでは、このラベルは0 (Z_¥rho = 0
  • そして、ラベル値Zは酔歩でランダムに決めるが、「どれくらい離れているときに、どれくらい遠くまで歩いているか」という意味合いで酔歩を決めるときE¥[(Z_a-Z_b)^2¥] = d_e(a,b)と書く。これは、T_e上の2点a,bについて、その間の「木的な距離d_e(a,b)」を酔歩する時間とみなしたときに、到達する座標変化がZ_a-Z_bの値になるように酔歩してZを決めよ、と言う意味である
  • 具体的には、木T_eの節ごとに酔歩を定めて、それを分岐部でつなぐことで、¥rhoの座標を0としたときの各所の座標Zを定める、というもの
  • さて。T_eとその上の値ラベルZが与えられた
  • 2点a,b ¥in T_e間に、以下のような「擬距離」を定める
  • 定めるにあたって2ステップを踏む
  • 第1ステップ
    • D^¥circ (a,b) = Z_a = Z_b - 2 max(min_{c ¥in [a,b} Z_c, min_{c \in [b,a]} Z_c)]
      • 木をaからbへ周回するか、bからaへ周回することにするとして、その周回にあたっては、わざわざ遠回りせずに、行き先に近い方の木のサイドから出発し、到着したらそれでOK、ということにする
      • そのようにするとしたうえで、周回中に出会う、最小のZ値のうち、大きめの方を採用して、そこを基準水位として、2点のZ値から基準水位まで降りた高さの和を、2点間のD^¥circ(a,b)とする、というもの
      • Z_a=Z_bであって、その周回路にZ_a=Z_bより小さい値がなければD^¥circ (a,b)=0である、ということであり、そのような関係の2点は「同一視」する、ということ
  • 第2ステップ
    • 上記のように定めたD^¥circ(a,b)を用いて、「本当の2点間の擬距離D^*(a,b)を次のように定める
    • 定義式には複数の書き方がある
      • 一つ目は木上の点ではなく、それに対応するブラウン散歩の時刻を使って与えるものである(複数の時刻が同一の木上の点に対応することに注意)
        • D^*(s,t) = ¥inf ¥{¥sum_{i=1}^k D^¥circ (s_i,t_i); k ¥ge 1, s=s1,t=t_k,d_e(t_i,s_{i+1})¥}
        • 木上の点 a,bについて、対応するブラウン散歩時刻を適当に取り(複数とれるとしても、同じaに対応する2時刻s1,s2についてはd_e(s1,s2)=0なる関係にあるので、上の式条件で「バイパスして距離を通算できる」ので問題ない
      • D^*(a,b) _{a,b¥in T_e} = ¥inf¥{ ¥sum_{i=1}^{k-1} D^¥circ (a_i,a_{i+1}) , k¥ge 1, a=a_1,a_2,....,a_k=b¥}
        • ここでは、木上の始点と終点とは固定されているが、その途中で何個の点を介してもよく(kの値は何でもよい)、その上で、ひたすら、下限をもたらす木上の点列を探索せよ、と言っている
    • 後掲のRmdの中では、「バイパスしてよい2点間」の距離を0にした上で、通ってよい点すべてをノードとする完全グラフの、グラフ上最短経路を計算することで、D^*を計算している

[][][][][]Rmdファイル『私のためのBrownian map構成法』 17:47

---
title: "ブラウニアンマップの構成"
author: "ryamada"
date: "2017年8月26日"
output: 
  html_document:
    toc: true
    toc_depth: 6
    number_section: true
---

```{r setup, include=FALSE}
library(rgl)
library(knitr)
knitr::opts_chunk$set(echo = TRUE)
knit_hooks$set(rgl = hook_rgl)
library(e1071)
library(ape)
library(igraph)
library(gtools)
```

# ブラウニアンマップの構成手順

ブラウニアンマップはS2同相なランダムな曲面。

このブラウニアンマップが、以下に示すランダムな木グラフ作成とそのグラフのノードへのランダムな数値ラベル付与をしたものに定義づけた距離関係情報の極限であることが知られている。

## ランダムな木グラフの作成

0からスタートし0に戻る1次元ブラウン運動のうち、0より大の値のみを取るそれを、ブラウン散歩(Brownian excursion)と呼ぶ。

このブラウン散歩から木グラフを構成することができる(構成法は後述する)。

この木グラフが、ブラウニアンマップ構成の1段階目である。

## 木構造を周回コースとしてのグラフ

作成したランダムな木グラフに沿って、ブラウン散歩を実行すると、すべてのエッジは両方向に1回ずつ、計2回、歩まれる。

エッジを2回歩くことを、エッジの両側面を歩くことだとみなすと、ランダム木グラフ上のブラウン散歩は、ルートノードから、木の周囲をぐるりと歩くことに相当する。

この周回コースは有向グラフとみなせる。

## ランダム数値のラベル付与

木グラフの直線成分ごとに酔歩をすることで、木グラフのノードにランダムな値が付与できる。

## 周回グラフ上の点の間の距離関係

周回グラフでは木グラフ上は同一の点であっても周回路の上では異なる点に相当することがある。
ブラウニアンマップの構成は、この周回グラフ上の点の間の距離関係を定めることにより実現される。

木グラフのノードに付与された値(この値は、周回グラフでは区別されていても、木の上で同一であれば、同じ値を持つことになる)と、木の構造とからノード間距離を定めるルールを導入する。このルールを定めることで、その木グラフが配置された二次元面に距離空間ができる。

この距離空間がブラウニアンマップである。

# ブラウン散歩によるランダムな木グラフの構成

## ブラウン散歩

0を出発して0に戻るブラウン運動はブラウニアンブリッジ(ブラウン橋)と言う。

ブラウン散歩のシミュレーション作成は、ブラウン橋を作成し、正の領域のみを通過するものができるまで、作成を繰り返すことで実現する。

```{r}
# Wiener bridge作成関数を持つパッケージ
#library(e1071)
my.rexcursion <- function(frequency = 1000, end = 1){
	succeeded <- FALSE
	while(!succeeded)
	{
		bridge <- rbridge(end = end, frequency = frequency)
		succeeded=all(bridge>=0)||all(bridge<=0)
	}
	return(c(0,c(abs(bridge))))
}
my.rwiener <- function(frequency = 1000, end = 1){
	c(0,cumsum(rnorm(end * frequency)/sqrt(frequency)))
}

fr <- 50
br.ex <- my.rexcursion(frequency=fr)
plot(br.ex,pch=20,type="b")
```

## ブラウン散歩から木を構成する

### ブラウン散歩上の点間距離の定義
ブラウン散歩で座標が大きくなるときには、未踏のエッジを作成し、
座標が小さくなるときには、来た道を戻ることを繰り返す。

このような木を作成するにあたり、以下の方法を採用する。

ブラウン散歩$e$上の2点$s,t$間の距離$d_e(s,t)$を以下のように定める。

$$
d_e(s,t) = x(s) + x(t) - 2 \times min_{u \in [s,t]} x(u)
$$

ただし、$x(u)$はパラメタ$u$でのブラウン散歩の座標とする。

ブラウン散歩上の点間の距離$d_e(s,t)$が定まった。
```{r}
my.Rtree.dist <- function(g,s,t){
	g[s] + g[t] - 2 * min(g[s:t])
}

my.Rtree.dist.mat <- function(g){
	L <- length(g)
	D.mat <- matrix(0,L,L)
	
	for(i in 1:(L-1)){
		for(j in (i+1):L){
			D.mat[i,j] <- D.mat[j,i] <- my.Rtree.dist(g,i,j)
		}
	}
	return(D.mat)
}
D.mat <- my.Rtree.dist.mat(br.ex)
image(D.mat)
```

### NJ法による木構造の作成
これを満足する木構造が存在することが知られているので、それを満足する木グラフをNJ法にて構成することができる。

```{r}
# nj法の関数を持つパッケージ
#library(ape)
tr <- nj(D.mat)

plot(nj(my.Rtree.dist.mat(br.ex)),type="u",show.tip.label=FALSE)
```

この木構造では、ブラウン散歩上の点に多数のノードが追加されるが、それらは、散歩上の点であることも知られているから、以下のような要領で、追加ノードを散歩上の点に対応するノードに変換する。

```{r}
# library(igraph)
my.extree <- function(br.ex,minval = 10^(-15)){
	n <- length(br.ex)
	D.mat <- my.Rtree.dist.mat(br.ex)
	tr <- nj(D.mat)
	el <- tr$edge
	el[which(el==n)] <- 1
	L <- tr$edge.length
	shortL <- (L < minval)
	while(max(el) > n){
		el <- t(apply(el,1,sort))
		tmp <- el - (n+0.5)
		tmp2 <- tmp[,1] * tmp[,2]
		tmp3 <- which(shortL & (tmp2 < 0))
		for(i in 1:length(tmp3)){
			el[which(el==el[tmp3[i],2])] <- el[tmp3[i],1]
		}
	}
	el.true <- which(el[,1]!=el[,2])
	el2 <- el[el.true,]
	L2 <- L[el.true]
	return(list(el=el2,len=L2))
}


extree <- my.extree(br.ex)
g <- graph.edgelist(extree$el,directed=FALSE)
g$weight <- extree$len
plot(g)
```


# 木の周り周回コースとしてのグラフ



```{r}
my.circuit.extree <- function(br.ex){
	extree <- my.extree(br.ex)
	g <- graph.edgelist(extree$el,directed=FALSE)
	n <- length(br.ex)
	ret <- c(1)
	for(i in 1:(n-2)){
		tmp <- shortest_paths(g,from=i,to=i+1)[[1]][[1]]
		ret <- c(ret,tmp[-1])
	}
	tmp <- shortest_paths(g,from=n-1,to=1)[[1]][[1]]
	ret <- c(ret,tmp[-1])
	ret[length(ret)] <- n
	return(ret)
}

circuit <- my.circuit.extree(br.ex)
plot(br.ex[circuit])
g.circuit <- graph.edgelist(cbind(circuit[1:(length(circuit)-1)],circuit[2:length(circuit)]))
plot(g.circuit) # 第1ノードと最終ノードは別ノードとして扱っているが、同一地点に相当することに注意
circuit
```


# ランダムな木グラフのノードへのランダムな値付与

木グラフのエッジの長さに応じて1次元ブラウン運動をさせ、それにより、エッジごとにエッジ始点からエッジ終点までの値の相対的増減値を発生させる。

その値を用いて、ルートノードの値を0としたときの各ノードの値を算出する。

```{r,rgl=TRUE}
my.edgeWiener <- function(lens,frequency=1000){
	ret <- rep(0,length(lens))
	for(i in 1:length(ret)){
		tmp <- rwiener(lens[i],frequency=frequency/lens[i])
		ret[i] <- tmp[length(tmp)]
	}
	return(ret)
}
my.edgeWiener(extree$len)

my.extreeZ <- function(br.ex){
	extree <- my.extree(br.ex)
	g <- graph.edgelist(extree$el,directed=FALSE)
	edgeZlen <- my.edgeWiener(extree$len)
	path.from.1 <- shortest_paths(g,from=1,output="epath")
	Z <- rep(0,length(br.ex)-1)
	for(i in 2:length(Z)){
		tmp <- as_ids(path.from.1[[2]][[i]])
		Z[i] <- sum(edgeZlen[tmp])
	}
	Z <- c(Z,0)
	return(Z)
}

Z <- my.extreeZ(br.ex)
plot3d(1:(fr+1),br.ex,Z,type="l")
```

# 周回グラフ上の点間距離

周回グラフ上の点間距離(擬距離)2段階の定義になっている。

## $D^{\circ}$

これは、ブラウン散歩のときに散歩座標から木グラフを作るときに導入した距離に似ている。

違いは、周回グラフであるので、順方向・逆方向の両方向を考慮して、短くなる方を採用するようになっていることである。


今、周回グラフがあって、ノード$p,q$に値$Z(p),Z(q)$がラベル付けされているとする。

この2点間に以下のような距離を定める。

$$
D^{\circ}(p,q) = Z(p) + Z(q) - 2 \times max(min_{r \in \{p \rightarrow q\}}Z(r),min_{r \in \{q \rightarrow p\}}Z(r))
$$
ただし$r\in \{p \rightarrow q\}$は$p$から$q$へと時計回りに巡回したときの周回グラフの部分を表し、
$r\in \{q \rightarrow p\}$は$q$から$p$へと時計回りに巡回したときの周回グラフの部分を表す。

2つの周回路での$Z$値の最小値のうち大きい方を取ることで、2点間の距離としては、「短め」になるような定義になっている。

$Z(p)=Z(q)$であって、どちらかの周回路での最小値が$Z(p)=Z(q)$であるときに、$D^{\circ}(p,q)=0$となることもわかる。

それ以外の場合は正の値を取る。

以下で、すべてのノード間の$D^{\circ}$を計算し、その行列を作成する。


```{r}
my.D.circuit <- function(v,Z,i,j){
	n <- length(v)
	Z. <- Z[v]
	vij <- i:j
	vji <- c(1:n,1:n)[j:(length(v)+i)]
	tmp <- max(min(Z.[vij]),min(c(Z.,Z.)[vji]))
	return(Z.[i]+Z.[j]-2*tmp)
}
# すべての
my.D.circuit.mat <- function(v,Z){
	n <- length(v)
	ret <- matrix(0,n,n)
	for(i in 1:(n-1)){
		for(j in (i+1):n){
			ret[i,j] <- ret[j,i] <- my.D.circuit(v,Z,i,j)
		}
	}
	return(ret)
}

Dcmat <- my.D.circuit.mat(circuit,Z)
image(Dcmat)
```

## $D(p,q)$

$D^{\circ}$ 自体を用いると、


$max(min_{r \in \{p \rightarrow q\}}Z(r),min_{r \in \{q \rightarrow p\}}Z(r))$の値が小さくなることがあり(周回路上に小さなZ値が出てくるため)、結果として、ある2点間の距離(擬距離)が長くなり過ぎて、三角不等式を満足することがなくなるなど、不具合があるという。

それを回避するために、$D^{\circ}$を使いつつ、異なる定義をする。

ブラウン散歩にて同一視された点s,tは$d_e(s,t)=0$の関係に関係にある、と表せるが、そのような点なら、バイパスしてもよいと言うルールで、以下のように定める。

$$
D(s,t) = Inf(\sum_{i=1}^k D^{\circ}(s_i,t_i)); s_1=s,t_k=k,d_e(s_{i+1},t_i)=0
$$

ただし、$k$の値は1でも良いし、他のどんな正の整数でもよいものとする。

$k=1$のときには、バイパスはせずに、$D^{\circ}(s,t)$を用いるし、$k>1$のときには、どこかしらの$d_e(s_{i+1},t_i)=0$な関係にある周回2点$t_i,s_{i+1}$でバイパスをする。

この$Inf$を取るにあたり、$s,t$の2点と、$d_e(s_{i+1},t_i)=0$を満足することから、バイパス候補になる点のみを取り出して、そのペア間に$D^{\circ}$の重みを持つ完全グラフを作成した上で、$s,t$間の最短距離を求める。

$d_e(s,t)$なる点の組み合わせを列挙する。

```{r}

Dcmat <- my.D.circuit.mat(circuit,Z)

Z <- my.extreeZ(br.ex)

```
まず、ブラウン散歩で下りから上りに転ずる点を識別する
```{r}
# d_e(s,t)=0のIDを取り出す
my.exUpDown <- function(br.ex){
	L <- length(br.ex)
	d <- diff(br.ex)
	#print(d)
	d.sign <- sign(d)
	#print(d.sign)
	d.d.sign <- diff(d.sign)
	dd.positive <- which(d.d.sign>0) + 1
	dd.negative <- which(d.d.sign<0) + 1
	#print(dd.positive)
	#print(dd.negative)
	col <- rep(1,L)
	col[dd.positive] <- 2
	col[dd.negative] <- 3
	#plot(br.ex,col=col,pch=20,type="b")
	return(col)
}
col <- my.exUpDown(br.ex)
plot(br.ex,col=col,pch=20,type="b")
plot(c(Z[circuit],Z[circuit]),type="b",pch=20,col=c(col[circuit],col[circuit]))
```

この識別情報から、周回路のノードIDで同等扱いするノードペアを列挙する
```{r}
valleyID <- which(my.exUpDown(br.ex)==2)
# d_e(s,t)に相当するcircuitノードペアを作る
my.valley.ids <- function(br.ex){
  valleyID <- which(my.exUpDown(br.ex)==2)
  el <- matrix(0,0,2)
  for(i in 1:length(valleyID)){
	  tmp <- which(circuit==valleyID[i])
	  tmp.comb <- combinations(length(tmp),2)
	  el <- rbind(el,cbind(tmp[tmp.comb[,1]],tmp[tmp.comb[,2]]))
  }
  return(list(id=unique(c(el)),pairs=el))
}
my.valley.ids(br.ex)
```

得られた周回路ノードと、着目2ノードとからなるサブグラフ(完全グラフ)を作り、最短距離を求める。

```{r}

#Z <- my.extreeZ(br.ex)

my.Dc.circuit <- function(i,j,Dcmat,valleys){
  tmp <- Dcmat[c(i,j,valleys$id),c(i,j,valleys$id)]
  n <- length(tmp[,1])
  pairs <- as.matrix(expand.grid(1:n,1:n))
  g <- graph.edgelist(pairs,directed=FALSE)
  ret <- distances(g,1,2,weights=c(tmp))
  return(ret)
}
my.Dc.circuit.mat <- function(br.ex,Z){
  valleys <- my.valley.ids(br.ex)
  circuit <- my.circuit.extree(br.ex)
  n <- length(circuit)
  ret <- matrix(0,n,n)
  Dcmat <- my.D.circuit.mat(circuit,Z)
  Dcmat2 <- Dcmat
  Dcmat2[(valleys$pairs[,1]-1) * n + valleys$pairs[,2]] <- 0

  for(i in 1:(n-1)){
    for(j in (i+1):n){
      Dcmat3 <- Dcmat2
      Dcmat3[i,] <- Dcmat[i,]
      Dcmat3[j,] <- Dcmat[j,]
      Dcmat3[,i] <- Dcmat[,i]
      Dcmat3[,j] <- Dcmat[,j]
      ret[i,j] <- ret[j,i] <- my.Dc.circuit(i,j,Dcmat3,valleys)
    }
  }
  ret
}
D <- my.Dc.circuit.mat(br.ex,Z)
image(D)
image(D==0)
```

# 検算

Z値が一番小さいノードからの距離は、単純に、Z値の差の絶対値になっている
```{r}
min.Z.id <- which(Z==min(Z))
circuit.min.Z.id <- which(circuit==min.Z.id)
for(i in 1:length(circuit.min.Z.id)){
  print(range(abs(Z[circuit]-Z[min.Z.id]) - D[circuit.min.Z.id[i],]))
  print(range(abs(Z[circuit]-Z[min.Z.id]) - Dcmat[circuit.min.Z.id[i],]))
}
```

周回路上、隣り合うノードの距離は、単純にZ値の差の絶対値
```{r}
for(i in 2:length(circuit)){
  print(D[i-1,i] - abs(Z[circuit[i-1]]-Z[circuit[i]]))
}
```

第一ノードと周回路最終ノードとの距離は0
```{r}
print(D[1,length(circuit)])
```

距離0のペアを検証する

```{r}
zero.pairs <- which(D==0,arr.ind=TRUE)
# 自身とのペアを除く
n <- length(circuit)
tmp <- zero.pairs[,1] != zero.pairs[,2]
zero.pairs. <- zero.pairs[tmp,]
zero.pairs. <- t(apply(zero.pairs.,1,sort))
for(i in 1:length(zero.pairs.[,1])){
  vij <- zero.pairs.[i,1]:zero.pairs.[i,2]
  vji <- c(1:n,1:n)[zero.pairs.[i,2]:(length(circuit)+zero.pairs.[i,1])]
  print(zero.pairs.[i,])
  print(circuit[zero.pairs.[i,]])
  print(min(Z[circuit[zero.pairs.[i,]]])-max(min(Z[circuit[vij]]),min(Z[circuit[vji]])))
}
```


$d_e(s1,s2)=d_e(s2,s3)=d_e(s3,s1)=0$なる3ノードが周回グラフには生じるが、すべてで$D(si,sj)=0$になるわけではないことが確認できる。
それと、「自身より大きなZ値なら迂回路を引ける」ということが対応することも確認できる。


トラックバック - http://d.hatena.ne.jp/ryamada/20170826