遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

BDDとZDD

BDD ZDD R
  • xは2^k行の真偽値表(２列)

library(igraph)
bzdd <- function(x,type="zdd"){
	v <- c(x)
	K <- log(length(v),2)
	# グラフに描くときに、判断分岐の高さごとにy軸座標を与えたいので、その値
	lay.y <- c()
	for(i in 1:K){
		lay.y <- c(lay.y,rep(i,2^(i-1)))
	}
	# 一番最後の群は「０，１」の真偽値の２ノード、その分のy軸座標
	lay.y <- c(lay.y,rep(K+1,2))
	g.tree <- graph.tree(sum(2^(0:(K-1))))
	# ２分岐木の隣接行列
	m.a <- as.matrix(get.adjacency(g.tree))
	# ただし、BDDでは、１ノードから出る２本のエッジに０，１のタイプ分けがある
	# それを隣接行列の値として1,2と区別する
	for(i in 1:length(m.a[,1])){
		m.a[i,which(m.a[i,]==1)[2]] <- 2
	}
	# 真偽値２ノードを付け加える
	m.a.2 <- matrix(0,length(m.a[,1])+2,length(m.a[1,])+2)
	m.a.2[1:length(m.a[,1]),1:length(m.a[1,])] <- m.a
	tmp.m <- matrix(v,ncol=2)
	# この真偽値情報を基にして、真偽値ノードへのエッジを隣接行列に加える
	# 最終判断項目が０でも１でもある真偽値になるなら、それは「１＋２＝３」の値を隣接行列に持たせる
	m.a.2[sum(2^(0:(K-2)))+which(tmp.m[,1]==0),length(m.a[1,])+1] <- 1
	m.a.2[sum(2^(0:(K-2)))+which(tmp.m[,2]==0),length(m.a[1,])+1] <- m.a.2[sum(2^(0:(K-2)))+which(tmp.m[,2]==0),length(m.a[1,])+1] + 2
	m.a.2[sum(2^(0:(K-2)))+which(tmp.m[,1]==1),length(m.a[1,])+2] <- 1
	m.a.2[sum(2^(0:(K-2)))+which(tmp.m[,2]==1),length(m.a[1,])+2] <- m.a.2[sum(2^(0:(K-2)))+which(tmp.m[,2]==1),length(m.a[1,])+2] + 2
	# 0(エッジなし)、１ Noのエッジ、2 Yesのエッジ、3 NoとYesの両方がそろったときのエッジ
	# BBD/ZDD化の本番
	#print(apply(m.a.2,1,sum))
	# ルールに従って隣接行列を更新する
	if(type=="zdd"){
		m.a.zdd <- m.a.2
		m.a.2 <- NULL
		#for(i in (K-1):1){
		for(i in (K-1):0){
			# 世代ごとに処理する。下流の世代から処理する
			# 今回の処理対象ノードを選ぶ
			if(i==0){
				ns <- c(1)
			}else{
				ns <- (sum(2^(0:(i-1))) + 1):sum(2^(0:i))
			}
			# この世代のノードを１個ずつ処理
			for(j in 1:length(ns)){
				# 処理中のノードのIDを取り出す
				current <- ns[j]
				tozero <- m.a.zdd[current,length(m.a.zdd[,1])-1]
				#print("tozero")
				#print(tozero)
				if(tozero == 2 | tozero == 3){
					for(jjj in 1:3){
						coming <- which(m.a.zdd[,current]==jjj)
						togo <- which(m.a.zdd[current,] == 1 | m.a.zdd[current,] == 3)
						#print("coming")
						#print(coming)
						#print("togo")
						#print(togo)
						if(length(coming)>0){
							m.a.zdd[coming,togo] <- m.a.zdd[coming,togo] + jjj
							m.a.zdd[coming,current] <- 0
							m.a.zdd[current,] <- 0
						}
						if(current == 1){
							m.a.zdd[current,]<- 0
						}
					}
				}
				#print(m.a.zdd)
			}
		}

	}else if(type=="bdd"){
		for(i in (K-1):1){
		#for(i in (K-1):0){
			# 世代ごとに処理する。下流の世代から処理する
			# 今回の処理対象ノードを選ぶ
			#if(i==0){
			#	ns <- c(1)
			#}else{
				ns <- (sum(2^(0:(i-1))) + 1):sum(2^(0:i))
			#}
			# この世代のノードを１個ずつ処理
			for(j in 1:length(ns)){
				# 処理中のノードのIDを取り出す
				current <- ns[j]
				#print(ns)
				# エッジが１本しか出ていなければ、そのノードは不要
				# 入ってくるエッジを出ていく先のノードに付け替える
				if(length(which(m.a.2[current,]>0))==1){
					togo <- which(m.a.2[current,]>0)
					# エッジの付け替えでは、エッジのタイプ値は変えない
					for(jj in 1:3){
						coming <- which(m.a.2[,current]==jj)
						m.a.2[coming,togo] <- m.a.2[coming,togo] + jj
						m.a.2[coming,current] <- 0
						m.a.2[current,] <- 0
					}
		# エッジが２本あるときは、同世代のノードに自分と同じパターンで下流につながっているならば、そのノードがあれば済むので、そのノードに付け替える
		# エッジのタイプ値については、上と同様に。
				}else if(length(which(m.a.2[current,]>1))==2){
					if(j < length(ns)-1){
						for(jj in (j+1):length(ns)){
							if(sum(abs(m.a.2[current,]-m.a.2[ns[jj],]))==0){
								togo <- ns[jj]
								for(jjj in 1:3){
									coming <- which(m.a.2[,current]==jjj)
									m.a.2[coming,togo] <- m.a.2[coming,togo] + jjj
									m.a.2[coming,current] <- 0
									m.a.2[current,] <- 0
								}
							}
						}
					}
					
				}
			}
		}

	}else if(type == "both"){
		m.a.zdd <- m.a.2
		#for(i in (K-1):1){
		for(i in (K-1):0){
			# 世代ごとに処理する。下流の世代から処理する
			# 今回の処理対象ノードを選ぶ
			if(i==0){
				ns <- c(1)
			}else{
				ns <- (sum(2^(0:(i-1))) + 1):sum(2^(0:i))
			}
			#ns <- (sum(2^(0:(i-1))) + 1):sum(2^(0:i))
			# この世代のノードを１個ずつ処理

			for(j in 1:length(ns)){
				# 処理中のノードのIDを取り出す
				current <- ns[j]
				#print(ns)
				if(length(ns)>1){
					# エッジが１本しか出ていなければ、そのノードは不要
					# 入ってくるエッジを出ていく先のノードに付け替える
					if(length(which(m.a.2[current,]>0))==1){
						togo <- which(m.a.2[current,]>0)
						# エッジの付け替えでは、エッジのタイプ値は変えない
						for(jj in 1:3){
							coming <- which(m.a.2[,current]==jj)
							m.a.2[coming,togo] <- m.a.2[coming,togo] + jj
							m.a.2[coming,current] <- 0
							m.a.2[current,] <- 0
						}
			# エッジが２本あるときは、同世代のノードに自分と同じパターンで下流につながっているならば、そのノードがあれば済むので、そのノードに付け替える
			# エッジのタイプ値については、上と同様に。
					}else if(length(which(m.a.2[current,]>1))==2){
						if(j < length(ns)-1){
							for(jj in (j+1):length(ns)){
								if(sum(abs(m.a.2[current,]-m.a.2[ns[jj],]))==0){
									togo <- ns[jj]
									for(jjj in 1:3){
										coming <- which(m.a.2[,current]==jjj)
										m.a.2[coming,togo] <- m.a.2[coming,togo] + jjj
										m.a.2[coming,current] <- 0
										m.a.2[current,] <- 0
									}
								}
							}
						}
						
					}

				}
				#print("current")
				#print(current)
				tozero <- m.a.zdd[current,length(m.a.zdd[,1])-1]
				#print("tozero")
				#print(tozero)
				if(tozero == 2 | tozero == 3){
					for(jjj in 1:3){
						coming <- which(m.a.zdd[,current]==jjj)
						togo <- which(m.a.zdd[current,] == 1 | m.a.zdd[current,] == 3)
						#print("coming")
						#print(coming)
						#print("togo")
						#print(togo)
						if(length(coming)>0){
							m.a.zdd[coming,togo] <- m.a.zdd[coming,togo] + jjj
							m.a.zdd[coming,current] <- 0
							m.a.zdd[current,] <- 0
						}
						if(current == 1){
							m.a.zdd[current,]<- 0
						}
					}
				}
				#print(m.a.zdd)
			}
		}

	}
	return(list(bdd.m = m.a.2,zdd.m = m.a.zdd))

}

x <- matrix(c(0,0,1,0,0,1,0,0),byrow=TRUE,ncol=2)
bzdd.out <- bzdd(x,type="both")
make.complex <- function(x){
	n <- length(x[,1])
	if(n==1){
		return(NULL)
	}
	m <- matrix(0,n,n)
	selected <- rep(1,n)
	for(i in 1:n){
		tmp <- apply(x,1,"-",x[i,])
		same <- apply((tmp==0),2,prod)
		bigger <- apply((tmp<=0),2,prod)
		smaller <- apply((tmp>=0),2,prod)
		if(length(which((smaller==1) & !(same==1)))>0){
			selected[i] <- 0
		}
		if(length(which(same[1:i]==1))>1){
			selected[i] <-0
		}
		m[i,which(bigger==1)] <- 2
		m[i,which(smaller==1)] <- 3
		m[i,which(same==1)] <- 1

	}
	tmp <- which(selected==1)
	if(length(tmp)==1){
		x <- matrix(x[which(selected==1),],nrow=1)
	}else{
		x <- x[which(selected==1),]
	}
	return(x)
}
n <- 8
Omega <- 1:n
k<- 5

x <- matrix(0,k,n)
for(i in 1:k){
	x[i,sample(Omega,sample(1:(n/2),1))] <- 1
}

x.c <- make.complex(x)


b <- as.matrix(expand.grid(rep(list(0:1),n)))
b <- b[,n:1]
x.c <- x.c[,n:1]
vv <- rep(0,length(b[,1]))
for(i in 1:length(b[,1])){
	tmp <- t(x.c) - b[i,]
	judge <- FALSE
	for(j in 1:length(tmp[1,])){
		if(length(which(tmp[,j]==-1))==0){
			judge <- TRUE
		}
	}
	if(judge){
		vv[i] <- 1
	}
}
x <- matrix(c(0,0,1,0,0,1,0,0),byrow=TRUE,ncol=2)
bzdd.out <- bzdd(x,type="both")
bdd <- graph.adjacency(sign(bzdd.out$bdd.m))
zdd <- graph.adjacency(sign(bzdd.out$zdd.m))
par(mfcol=c(1,2))
plot(bdd)
plot(zdd)
par(mfcol=c(1,1))