グラフのゼータ関数 その後

• 資料
• グラフのvertex zeta 関数はノード数ｘノード数の行列を使って以下のように表せる
  • 
    • ただし、uは複素数、mとnはグラフGのエッジ数とノード数、はMの行列式で、Iはノード数ｘノード数の単位行列、AはGの隣接行列、DはGのノードの次数を対角成分とする対角行列
• また、同じvertex zeta 関数はエッジ数ｘエッジ数の行列Wを使って
  • 
    • ただし、uは複素数、はMの行列式で、Wは(エッジ数x2)ｘ(エッジ数x2)の行列である。グラフは無向グラフであり、そのエッジに2方向を考慮して2エッジを対象とする。今、あるエッジの終点があるエッジの始点であるとき、その2エッジのペアに対応するWのセルに１を立て、他は０とする。ただし、始点と終点を取り換えただけのエッジペアに対応するセルには０を立てる
  • さらに。エッジ長がすべて１のときは、上記でよいが、エッジ長が１とは限らないときには工夫が必要となる。エッジ長がすべて整数であるときには、その長さに応じてノードを追加し、すべてのエッジが１のグラフに変換すればよい
• ちなみに、エッジの行列の行列式でさらっとvertex zeta 関数になるのは：
  • エッジ連結について、全部１が立っているとパスはそれを掛け合わせると、相変わらず１で、それが「存在する」という情報になるから。したがって、そこにひとしなみに複素数uを掛けるとの項がエッジ数だけ冪乗されて現れる。これは、行列Wの要素をuとすることに同じ
  • さらに、エッジに長さ重みを入れたいときには、この原理を使って、行列要素をとエッジごとに変えることで実現できる。ちなみにはエッジペア(a,b)のそれぞれのエッジの重み
• さらに小さな行列の行列式計算に変形することもできる。それがpath zeta 関数
  • グラフに全域木を取り出し、全域木に含まれないエッジを取り出す。この全域木に含まれないエッジ数の2倍の二乗のサイズの行列の行列式計算にできる
  • 補全域木グラフの両方向を考えたエッジの集合を考える
  • それらのペアに対して、自身の逆方向にあたるエッジとのペアのときには、対応する行列成ー分を０にする
  • それ以外の場合には、ゼータ関数のuに重みを考慮した値を与える
  • その値は、補全域木グラフの２つのエッジが接続した木の2点間について、木の上のパス(最短)について、最短距離を取る。edge zeta 関数ではつながるエッジペアについて、２エッジの長さの平均を用いた。それと同様にエッジのつながりについて足し合わせて行きたいが、エッジペアをタンデムに足し合わせていくことは、結局そのグラフ距離を計算することのようなもの
  • こちらの９章にmathematica のコードｇあるのでそれを読むのが早そう

• 以下は、ノード数ｘノード数の行列でのゼータ関数、(エッジ数x2)ｘ(エッジ数x2)の行列でのゼータ関数、エッジ長に合わせてノードを補う関数である

library(igraph)
library(complexplus)
my.Ihara.zeta.elem <- function(g){
	A <- as.matrix(as_adjacency_matrix(g))
	n.v <- length(V(g))
	n.e <- length(E(g))
	r <- n.e - n.v
	D <- diag(degree(g))
	H <- diag(degree(g)-1) 
	I <- diag(rep(1,n.v))
	return(list(r=r,A=A,H=H,I=I,D=D,n.v=n.v,n.e=n.e))
}
my.Ihara.zeta.ori <- function(g,u){
	elem <- my.Ihara.zeta.elem(g)
	1/((1-u^2)^elem$r * Det(elem$I-u*elem$A + u^2*elem$H))
}

my.Bartholdi.zeta <- function(g,u,t){
	elem <- my.Ihara.zeta.elem(g)
	tmp <- (1-(1-u)^2*t^2)^elem$r * Det(elem$I-t*elem$A+(1-u)*(elem$D-(1-u)*elem$I)*t^2)
	return(1/tmp)
}

my.Ihara.zeta <- function(g,u){
	my.Bartholdi.zeta(g,0,u)
}
# (エッジ長ｘ２）^2 行列でのゼータ関数
my.Ihara.zeta.e <- function(g,u){
	el <- as_edgelist(g)
	el2 <- rbind(el,cbind(el[,2],el[,1]))
	n.e <- length(el[,1])
	edge.mat <- matrix(0,n.e*2,n.e*2)
	for(i in 1:(n.e*2)){
		st <- el2[i,1]
		ed <- el2[i,2]
		tmp <- which(el2[,1]==ed & el2[,2]!=st)
		edge.mat[i,tmp] <- 1
	}
	I <- diag(rep(1,n.e*2))
	tmpmat <- I - edge.mat * u
	return(1/Det(tmpmat))
}

# エッジに整数長さがあるときに、すべてのエッジ長が１となるように
# ノードを加えたグラフに変換する
# E(g)$weight <- sample(1:20,n.e)

my.graph.inflation.weight <- function(g){
	w <- E(g)$weight
	el <- as_edgelist(g)
	A <- as.matrix(as_adjacency_matrix(g))
	for(i in 1:length(w)){
		if(w[i] > 1){
			nv <- w[i]-1
			st <- el[i,1]
			ed <- el[i,2]
			A[st,ed] <- A[ed,st] <- 0
			n <- length(A[,1])
			newA <- matrix(0,n+nv,n+nv)
			newA[1:n,1:n] <- A
			newA[n+1,st] <- newA[st,n+1] <- 1
			newA[n+nv,ed] <- newA[ed,n+nv] <- 1
			if(nv > 1){
				for(j in 1:(nv-1)){
					newA[n+j,n+j+1] <- newA[n+j+1,n+j] <-1
				}
			}
			A <- newA
		}
	}
	newg <- graph.adjacency(A,mode="undirected")
	return(newg)
}

E(g)$weight <- rep(1,n.e)
newg <- my.graph.inflation.weight(g)


fu <- function(u){
	tmp <- (-1)*(3*u-1)*(u+1)^5*(u-1)^6*(3*u^2+u+1)^4
	return(1/tmp)
}

#g <- sample_gnp(10, 2/10)
n <- 5
adj <- matrix(1,n,n)
diag(adj) <- 0
g <- graph.adjacency(adj,mode="undirected")

x <- y <- seq(from=-0.8,to=0.8,length=100)

xy <- expand.grid(x,y)
xy. <- xy[,1] + 1i*xy[,2]

vs <- rep(0,length(xy.))
vs2 <- vs
for(i in 1:length(vs)){
	vs[i] <- my.Ihara.zeta(g,xy.[i])
	vs2[i] <- fu(xy.[i])
}
my.Ihara.zeta(g,vs[55])
fu(vs[55])
plot(xy,col=abs(Mod(vs)))


fu2 <- function(u){
	tmp <- (1-u^10)^5*(1-3*u^5)*(1-u^5)*(1+u^5+3*u^10)
	return(1/tmp)
}


fu(vs[515]^5)
fu2(vs[515])