空間が広い - ryamadaの遺伝学・遺伝統計学メモ

こちらで数学知らずの乳幼児による度数分布平滑化手法というのをやってみた
多次元空間の場合も、直交軸に関する２階の差分を評価することで、いい感じの平滑化を取り出すことができるようだった
今度は、点の数がものすごく多いとき、どうするかと言う話
生物のように、発熱もしないエネルギー効率がものすごく優秀な光学センサーと情報処理系を持っていて、超高度な並列処理回路になっていれば問題ないが、お手軽PCでそれをまねようとすると、空間の広さ(度数をとる単位超立方体の数)が問題になる
度数分布を考えるとき、すべての超立方体についてそこに存在する点の数を知ろうとするのが普通だろうが、空間が広くなると、結構粗い区画(１軸50分割とか)でも、超立方体の数は $50^n$ になり、nが10とかになるとやりきれない数になる。 $50^10$ に比べれば、大量のシグナル( $10^5$ など)はかわいいもので、シグナル件数を相手にする方がまだましになってくる。それがたとえ、 $(10^5)^2$ というように、全シグナルの間のペアワイズな評価を必要とするとしても、同程度の厄介さかげんになってくる。空間の区画数を相手にするか、シグナル点を相手にするかは、場合によるが、どちらもできる方がよいだろう

# データの座標によらず、視野内に収めるための補正
vision.scale <- function(x,scl){
	(x-min(x))/(max(x)-min(x)) * scl + (1-scl)/2
}
# 観察対象をきちんと視野内に収めるために、周辺部にはシグナルがないように調整する
vision.standard <- function(x,scl = 0.9){
	rg <- apply(x,2,range)
	tmp <- apply(x,2,vision.scale, scl)
	return(list(st.x = tmp,rg=rg))
}
# 視野全体のピクセル数を指定して、該当する超立方体の番地にする
round.coords <- function(x,n.cell=50){
	d <- length(x[1,])
	matrix(floor(x * n.cell),ncol=d)
}
library(igraph)
cluster.g <- function(d){
	sign.d <- sign(-(sign(d)-1))
	g <- graph.adjacency(sign.d)
	clusters(g)
}
unique.dist <- function(d){
	cl <- cluster.g(d)
	unique.id <- which(!duplicated(cl[[1]]))
	return(list(cl=cl,unique.id=unique.id))
}
one.direction.2 <- function(unique.d,tmp.x){
	d <- length(tmp.x[1,])
	# 距離１のエッジ
	### １エッジが２回登場している…
	up.ud <- unique.d
	up.ud[upper.tri(up.ud)] <- 0
	ones <- matrix(which(up.ud==1,arr.ind=TRUE),ncol=2)
	# 距離２のエッジ
	twos <- matrix(which(up.ud==2,arr.ind=TRUE),ncol=2)
	# 距離１のエッジの軸と向きを確認
	# 距離２が同方向のそれであることを確認し、その軸と向きを確認
	#tmp.x <- round.x[[i]][unique.id,]
	diff.ones <- matrix(tmp.x[ones[,1],],ncol=d) - matrix(tmp.x[ones[,2],],ncol=d)
	ones.axis <- apply(abs(diff.ones),1,which.max)
	ones.dir <- apply(diff.ones,1,sum)
	# 軸ごとに線分状グラフ(複数)を作成
	ones.directed <- ones
	ones.directed[which(ones.dir==-1),1] <- ones[which(ones.dir==-1),2]
	ones.directed[which(ones.dir==-1),2] <- ones[which(ones.dir==-1),1]
	trees <- list()
	for(i in 1:d){
		tmp.dir <- which(ones.axis==i)
		if(length(tmp.dir)>0){
			trees[[i]] <- tandem.one(matrix(ones.directed[tmp.dir,],ncol=2))
		}
		
	}
	#trees <- tandem.one(ones.directed)
	
	diff.twos <- matrix(tmp.x[twos[,1],],ncol=d) - matrix(tmp.x[twos[,2],],ncol=d)
	sign.diff.twos <- sign(abs(diff.twos))
	twos.line.id <- which(apply(sign.diff.twos,1,sum)==1)
	if(length(twos.line.id)>0){
		twos.line <- twos[twos.line.id,]
		diff.twos.line <- diff.twos[twos.line.id,]
		twos.axis <- apply(abs(diff.twos.line),1,which.max)
		twos.dir <- sign(apply(diff.twos.line,1,sum))
		twos.line.directed <- twos.line
		twos.line.directed[which(twos.dir==-1),1] <- twos.line[which(twos.dir==-1),2]
		twos.line.directed[which(twos.dir==-1),2] <- twos.line[which(twos.dir==-1),1]
		#nec.unnec.two <- omit.twos(cbind(ones[,st.one],ones[,end.one]),cbind(twos.line[,st.two],twos.line[,end.two]))
		for(i in 1:d){
			tmp.dir <- which(twos.axis==i)
			if(length(tmp.dir)>0){
				trees[[i]] <- tandem.two(trees[[i]],matrix(twos.line.directed[tmp.dir,],ncol=2))
			}
		}

	}

	return(list(trees))
}
one.direction <- function(unique.id,unique.d,tmp.x){
	d <- length(tmp.x[1,])
	# 距離１のエッジ
	### １エッジが２回登場している…
	up.ud <- unique.d
	up.ud[upper.tri(up.ud)] <- 0
	ones <- matrix(which(up.ud==1,arr.ind=TRUE),ncol=2)
	# 距離２のエッジ
	twos <- matrix(which(up.ud==2,arr.ind=TRUE),ncol=2)
	# 距離１のエッジの軸と向きを確認
	# 距離２が同方向のそれであることを確認し、その軸と向きを確認
	#tmp.x <- round.x[[i]][unique.id,]
	diff.ones <- matrix(tmp.x[ones[,1],],ncol=d) - matrix(tmp.x[ones[,2],],ncol=d)
	ones.axis <- apply(abs(diff.ones),1,which.max)
	ones.dir <- apply(diff.ones,1,sum)
	# 軸ごとに線分状グラフ(複数)を作成
	ones.directed <- ones
	ones.directed[which(ones.dir==-1),1] <- ones[which(ones.dir==-1),2]
	ones.directed[which(ones.dir==-1),2] <- ones[which(ones.dir==-1),1]
	trees <- list()
	for(i in 1:d){
		tmp.dir <- which(ones.axis==i)
		if(length(tmp.dir)>0){
			trees[[i]] <- tandem.one(matrix(ones.directed[tmp.dir,],ncol=2))
		}
		
	}
	#trees <- tandem.one(ones.directed)
	
	diff.twos <- matrix(tmp.x[twos[,1],],ncol=d) - matrix(tmp.x[twos[,2],],ncol=d)
	sign.diff.twos <- sign(abs(diff.twos))
	twos.line.id <- which(apply(sign.diff.twos,1,sum)==1)
	if(length(twos.line.id)>0){
		twos.line <- twos[twos.line.id,]
		diff.twos.line <- diff.twos[twos.line.id,]
		twos.axis <- apply(abs(diff.twos.line),1,which.max)
		twos.dir <- sign(apply(diff.twos.line,1,sum))
		twos.line.directed <- twos.line
		twos.line.directed[which(twos.dir==-1),1] <- twos.line[which(twos.dir==-1),2]
		twos.line.directed[which(twos.dir==-1),2] <- twos.line[which(twos.dir==-1),1]
		#nec.unnec.two <- omit.twos(cbind(ones[,st.one],ones[,end.one]),cbind(twos.line[,st.two],twos.line[,end.two]))
		for(i in 1:d){
			tmp.dir <- which(twos.axis==i)
			if(length(tmp.dir)>0){
				trees[[i]] <- tandem.two(trees[[i]],matrix(twos.line.directed[tmp.dir,],ncol=2))
			}
		}

	}

	return(list(trees))
}
omit.twos <- function(ones,twos){
	nec <- unnec <- rep(0,length(twos[,1]))
	for(i in 1:length(twos[,1])){
		if(sum(which(ones[,1] == twos[i,1])) * sum(which(ones[,2] == twos[i,2]))){
			unnec[i] <- 1
		}else{
			nec[i] <- 1
		}
	}
	return(list(nec=nec,unnec=unnec))
}
tandem.one <- function(m){
	#print(m)
	trees <- list(m[1,])
	st.end <- matrix(m[1,],nrow=1)
	if(length(m[,1])>1){
		for(i in 2:length(m[,1])){
			#print(i)
			#print(st.end)
			tailer <- which(st.end[,2]==m[i,1])
			header <- which(st.end[,1]==m[i,2])
			#print("here")
			if(length(tailer)>0){
				if(length(header)>0){
					trees[[tailer]] <- c(trees[[tailer]],trees[[header]])
					trees[[header]] <- c()
					st.end[tailer,2] <- st.end[header,2]
					st.end <- matrix(st.end[-header,],ncol=2)
					#print(trees)
					#print(st.end)
				}else{
					trees[[tailer]] <- c(trees[[tailer]],m[i,2])
					st.end[tailer,2] <- m[i,2]
				}
			}else{
				if (length(header)>0){
					trees[[header]] <- c(m[i,1],trees[[header]])
					st.end[header,1] <- m[i,1]
				}else{
					trees[[length(st.end[,1])+1]] <- m[i,]
					st.end <- rbind(st.end,matrix(m[i,],nrow=1))
				}
			}
		}
	}
	return(list(trees=trees,st.end=st.end))
}
tandem.two <- function(trees,m){
	
	tr <- trees$trees
	st.end <- trees$st.end
	#print(trees)
	#print(tr)
	#print(st.end)
	for(i in 1:length(m[,1])){
		tailer <- which(st.end[,2]==m[i,1])
		header <- which(st.end[,1]==m[i,2])
		if(length(tailer)>0){
			if(length(header)>0){
				tr[[tailer]] <- c(tr[[tailer]],0,tr[[header]])
				tr[[header]] <- c()
				st.end[tailer,2] <- st.end[header,2]
				st.end <- matrix(st.end[-header,],ncol=2)
			}else{
				tr[[tailer]] <- c(tr[[tailer]],0,m[i,2])
				st.end[tailer,2] <- m[i,2]
			}
		}else{
			if(length(header)>0){
				tr[[header]] <- c(m[i,1],0,tr[[header]])
				st.end[header,1] <- m[i,1]
			}else{
				if(length(which(unlist(tr)==m[i,1]))>0 & length(which(unlist(tr)==m[i,2]))>0){
					
				}else{
					tr[[length(st.end[,1])+1]] <- c(m[i,1],0,m[i,2])
					st.end <- rbind(st.end,matrix(m[i,],nrow=1))
				}
				
			}
		}
	}
	return(list(trees=tr,st.end=st.end))
}


d <- 3
x1 <- c(rnorm(1000,10),rnorm(500,80,20),runif(1500,30,50))
x2 <- c(rnorm(500,5),rnorm(1500,20,10),runif(1000,4,50))
x3 <- c(rnorm(1500,40),rnorm(500,10,100),runif(1000,40,50))
x <- cbind(x1,x2,x3)
x.st <- vision.standard(x,scl=0.7)
x <- x.st$st.x
n <- length(x[,1])

library(rgl)
#plot3d(x)

Ls <- seq(from=min(dist(x,method="manhattan")),to=1,length=10)
Ls <- Ls[-1]
Ls <- Ls[-1]

round.x <- list()
dist.out <- dist.m <- cl <- uni.d <- trees <- list()
par(ask=FALSE)
for(i in 1:length(Ls)){
	round.x[[i]] <- round.coords(x,1/Ls[i])
	dist.out[[i]] <- dist(round.x[[i]],method="manhattan")
	#cl[[i]] <- cluster.g(dist.out[[i]])
	tmp <- unique.dist(dist.out[[i]])
	cl[[i]] <- tmp$cl
	unique.id <- tmp$unique.id
	unique.d <- as.matrix(dist.out[[i]])[unique.id,unique.id]
	tmp.x <- matrix(round.x[[i]][unique.id,],ncol=d)
	
	trees[[i]] <- one.direction(unique.id,unique.d,tmp.x)
	
	#plot3d(round.x[[i]])
	#image(dist.m[[i]])
	#plot(sort(dist.m[[i]]))
}

# データの座標によらず、視野内に収めるための補正
vision.scale <- function(x,scl){
	(x-min(x))/(max(x)-min(x)) * scl + (1-scl)/2
}
# 観察対象をきちんと視野内に収めるために、周辺部にはシグナルがないように調整する
vision.standard <- function(x,scl = 0.9){
	rg <- apply(x,2,range)
	tmp <- apply(x,2,vision.scale, scl)
	return(list(st.x = tmp,rg=rg))
}
# 視野全体のピクセル数を指定して、該当する超立方体の番地にする
round.coords <- function(x,n.cell=50){
	d <- length(x[1,])
	matrix(floor(x * n.cell),ncol=d)
}

d <- 3
x1 <- c(rnorm(1000,10),rnorm(500,80,20),runif(1500,30,50))
x2 <- c(rnorm(500,5),rnorm(1500,20,10),runif(1000,4,50))
x3 <- c(rnorm(1500,40),rnorm(500,10,100),runif(1000,40,50))
x <- cbind(x1,x2,x3)
x.st <- vision.standard(x,scl=0.7)
x <- x.st$st.x
n <- length(x[,1])

library(rgl)
#plot3d(x)

round.x <- count.x <- membership.x <- cluster.x <- represent.x <- infos <- list()
dist.mat.x <- list()
trees.out <- list()
dist.mat.x[[1]] <- NULL
round.x[[1]] <- x
N <- length(x[,1])
membership.x[[1]] <- 1:N
count.x[[1]] <- rep(1,N)
cluster.x[[1]] <- 1:N
represent.x[[1]] <- 1:N

Ls <- 1/(c(16,8,4,2))

for(i in 1:length(Ls)){
	if(length(matrix(round.x[[i]],ncol=d)[,1])>1){
		if(i==1){
			tmp <- round.coords(x,1/Ls[i])
		}else{
			tmp <- round.coords(round.x[[i]]*Ls[i-1],1/Ls[i])
		}
		
		n <- length(tmp[,1])
		membership.x[[i+1]] <- rep(0,n)
		count.x[[i+1]] <- vector()
		represent.x[[i+1]] <- vector()
		tobechecked <- rep(1,n)
		while(sum(tobechecked)>0){
			current.checker <- which(tobechecked==1)[1]
			current.tobechecked <- which(tobechecked==1)[-1]
			#if(length(current.tobechecked)>0){
				diffs <- matrix(t(tmp[current.tobechecked,])-tmp[current.checker,],nrow=d)
				diff.all <- apply(diffs,2,sum)
				zeros <- which(diff.all==0)
				membership.x[[i+1]][c(current.checker,current.tobechecked[zeros])] <- max(membership.x[[i+1]])+1
				
				tobechecked[c(current.checker,current.tobechecked[zeros])] <- 0
				represent.x[[i+1]] <- c(represent.x[[i+1]],current.checker)
				count.x[[i+1]] <- c(count.x[[i+1]],sum(count.x[[i]][c(current.checker,current.tobechecked[zeros])]))
		}
		round.x[[i+1]] <- tmp[represent.x[[i+1]],]
		dist.mat.x[[i+1]] <- as.matrix(dist(round.x[[i+1]],method="manhattan"))
		trees.out[[i+1]] <- one.direction.2(dist.mat.x[[i+1]],round.x[[i+1]])
		tmp.prob <- count.x[[i+1]]/sum(count.x[[i+1]])
		infos[[i+1]] <- sum(tmp.prob*log(tmp.prob))
	}else{
		break
	}
}

平滑化のためには、次元方向の２階の差分について考える必要がある
- 空間が広いので、ほとんどの超立方体では、その周辺に関して０の平原になっているとすれば、シグナルが帰属している超立方体とその隣接超立方体についてのみ２階差分を取ればよい
- しかも、今、次元の数の軸について、軸の組合せは気にせず、軸の方向のみを考えるとき、ある超立方体に関する２階の差分は、すべての隣の超立方体に帰属するシグナル数だけがわかればよい
- また、軸方向についてのみを考えるときには、軸方向ごとに処理を分離することができるのも、ハンドリング上ありがたい特徴となる
- シグナルが帰属する超立方体はすでに上記の処理で検出されている。このほかに必要なのは、その隣接超立方体であって、シグナルが帰属していないものである
- そのようなシグナルなし超立方体は２つに分けることができて、ある方向に関して両隣りにシグナルがあるか、型隣りにしかシグナルがないか、のいずれかである
- ここで、シグナルなし超立方体に関して全方向を同時に考慮するとなると、場合分けが非常に多くなるが、軸方向別に分離して処理できることから、ある超立方体を考えるときに、ある軸についての両隣りの状態のみを考慮して、それを軸数分だけ積み上げることで事足りる
- このように、シグナルなし超立方体を扱うためには、このような超立方体を２つのシグナルありの超立方体がはさんでいるかそうでないかの評価は必要となる。従って、シグナルあり超立方体同士の関係としては、隣(ハミング距離が１)か否かと、ある方向に距離２(折れ曲がっての距離ではなく)の飛び隣か否かの２点のみの確認が必要である
- その上で、軸別に、「隣」「飛び隣」の具合での直線上グラフを１個以上作成し、それについて２階差分を合算することとなる

# データの座標によらず、視野内に収めるための補正
vision.scale <- function(x,scl){
	(x-min(x))/(max(x)-min(x)) * scl + (1-scl)/2
}
# 観察対象をきちんと視野内に収めるために、周辺部にはシグナルがないように調整する
vision.standard <- function(x,scl = 0.9){
	rg <- apply(x,2,range)
	tmp <- apply(x,2,vision.scale, scl)
	return(list(st.x = tmp,rg=rg))
}
# 視野全体のピクセル数を指定して、該当する超立方体の番地にする
round.coords <- function(x,n.cell=50){
	d <- length(x[1,])
	matrix(floor(x * n.cell),ncol=d)
}
library(igraph)
cluster.g <- function(d){
	sign.d <- sign(-(sign(d)-1))
	g <- graph.adjacency(sign.d)
	clusters(g)
}
unique.dist <- function(d){
	cl <- cluster.g(d)
	unique.id <- which(!duplicated(cl[[1]]))
	return(list(cl=cl,unique.id=unique.id))
}

one.direction <- function(unique.id,unique.d,tmp.x){
	d <- length(tmp.x[1,])
	# 距離１のエッジ
	### １エッジが２回登場している…
	up.ud <- unique.d
	up.ud[upper.tri(up.ud)] <- 0
	ones <- matrix(which(up.ud==1,arr.ind=TRUE),ncol=2)
	# 距離２のエッジ
	twos <- matrix(which(up.ud==2,arr.ind=TRUE),ncol=2)
	# 距離１のエッジの軸と向きを確認
	# 距離２が同方向のそれであることを確認し、その軸と向きを確認
	#tmp.x <- round.x[[i]][unique.id,]
	diff.ones <- matrix(tmp.x[ones[,1],],ncol=d) - matrix(tmp.x[ones[,2],],ncol=d)
	ones.axis <- apply(abs(diff.ones),1,which.max)
	ones.dir <- apply(diff.ones,1,sum)
	# 軸ごとに線分状グラフ(複数)を作成
	ones.directed <- ones
	ones.directed[which(ones.dir==-1),1] <- ones[which(ones.dir==-1),2]
	ones.directed[which(ones.dir==-1),2] <- ones[which(ones.dir==-1),1]
	trees <- list()
	for(i in 1:d){
		tmp.dir <- which(ones.axis==i)
		if(length(tmp.dir)>0){
			trees[[i]] <- tandem.one(matrix(ones.directed[tmp.dir,],ncol=2))
		}
		
	}
	#trees <- tandem.one(ones.directed)
	
	diff.twos <- matrix(tmp.x[twos[,1],],ncol=d) - matrix(tmp.x[twos[,2],],ncol=d)
	sign.diff.twos <- sign(abs(diff.twos))
	twos.line.id <- which(apply(sign.diff.twos,1,sum)==1)
	if(length(twos.line.id)>0){
		twos.line <- twos[twos.line.id,]
		diff.twos.line <- diff.twos[twos.line.id,]
		twos.axis <- apply(abs(diff.twos.line),1,which.max)
		twos.dir <- sign(apply(diff.twos.line,1,sum))
		twos.line.directed <- twos.line
		twos.line.directed[which(twos.dir==-1),1] <- twos.line[which(twos.dir==-1),2]
		twos.line.directed[which(twos.dir==-1),2] <- twos.line[which(twos.dir==-1),1]
		#nec.unnec.two <- omit.twos(cbind(ones[,st.one],ones[,end.one]),cbind(twos.line[,st.two],twos.line[,end.two]))
		for(i in 1:d){
			tmp.dir <- which(twos.axis==i)
			if(length(tmp.dir)>0){
				trees[[i]] <- tandem.two(trees[[i]],matrix(twos.line.directed[tmp.dir,],ncol=2))
			}
		}

	}

	return(list(trees))
}
omit.twos <- function(ones,twos){
	nec <- unnec <- rep(0,length(twos[,1]))
	for(i in 1:length(twos[,1])){
		if(sum(which(ones[,1] == twos[i,1])) * sum(which(ones[,2] == twos[i,2]))){
			unnec[i] <- 1
		}else{
			nec[i] <- 1
		}
	}
	return(list(nec=nec,unnec=unnec))
}
tandem.one <- function(m){
	#print(m)
	trees <- list(m[1,])
	st.end <- matrix(m[1,],nrow=1)
	if(length(m[,1])>1){
		for(i in 2:length(m[,1])){
			#print(i)
			#print(st.end)
			tailer <- which(st.end[,2]==m[i,1])
			header <- which(st.end[,1]==m[i,2])
			#print("here")
			if(length(tailer)>0){
				if(length(header)>0){
					trees[[tailer]] <- c(trees[[tailer]],trees[[header]])
					trees[[header]] <- c()
					st.end[tailer,2] <- st.end[header,2]
					st.end <- matrix(st.end[-header,],ncol=2)
					#print(trees)
					#print(st.end)
				}else{
					trees[[tailer]] <- c(trees[[tailer]],m[i,2])
					st.end[tailer,2] <- m[i,2]
				}
			}else{
				if (length(header)>0){
					trees[[header]] <- c(m[i,1],trees[[header]])
					st.end[header,1] <- m[i,1]
				}else{
					trees[[length(st.end[,1])+1]] <- m[i,]
					st.end <- rbind(st.end,matrix(m[i,],nrow=1))
				}
			}
		}
	}
	return(list(trees=trees,st.end=st.end))
}
tandem.two <- function(trees,m){
	
	tr <- trees$trees
	st.end <- trees$st.end
	#print(trees)
	#print(tr)
	#print(st.end)
	for(i in 1:length(m[,1])){
		tailer <- which(st.end[,2]==m[i,1])
		header <- which(st.end[,1]==m[i,2])
		if(length(tailer)>0){
			if(length(header)>0){
				tr[[tailer]] <- c(tr[[tailer]],0,tr[[header]])
				tr[[header]] <- c()
				st.end[tailer,2] <- st.end[header,2]
				st.end <- matrix(st.end[-header,],ncol=2)
			}else{
				tr[[tailer]] <- c(tr[[tailer]],0,m[i,2])
				st.end[tailer,2] <- m[i,2]
			}
		}else{
			if(length(header)>0){
				tr[[header]] <- c(m[i,1],0,tr[[header]])
				st.end[header,1] <- m[i,1]
			}else{
				if(length(which(unlist(tr)==m[i,1]))>0 & length(which(unlist(tr)==m[i,2]))>0){
					
				}else{
					tr[[length(st.end[,1])+1]] <- c(m[i,1],0,m[i,2])
					st.end <- rbind(st.end,matrix(m[i,],nrow=1))
				}
				
			}
		}
	}
	return(list(trees=tr,st.end=st.end))
}


d <- 3
x1 <- c(rnorm(1000,10),rnorm(500,80,20),runif(1500,30,50))
x2 <- c(rnorm(500,5),rnorm(1500,20,10),runif(1000,4,50))
x3 <- c(rnorm(1500,40),rnorm(500,10,100),runif(1000,40,50))
x <- cbind(x1,x2,x3)
x.st <- vision.standard(x,scl=0.7)
x <- x.st$st.x
n <- length(x[,1])

library(rgl)
#plot3d(x)

Ls <- seq(from=min(dist(x,method="manhattan")),to=1,length=10)
Ls <- Ls[-1]
Ls <- Ls[-1]

round.x <- list()
dist.out <- dist.m <- cl <- uni.d <- trees <- list()
par(ask=FALSE)
for(i in 1:length(Ls)){
	round.x[[i]] <- round.coords(x,1/Ls[i])
	dist.out[[i]] <- dist(round.x[[i]],method="manhattan")
	#cl[[i]] <- cluster.g(dist.out[[i]])
	tmp <- unique.dist(dist.out[[i]])
	cl[[i]] <- tmp$cl
	unique.id <- tmp$unique.id
	unique.d <- as.matrix(dist.out[[i]])[unique.id,unique.id]
	tmp.x <- matrix(round.x[[i]][unique.id,],ncol=d)
	
	trees[[i]] <- one.direction(unique.id,unique.d,tmp.x)
	
	#plot3d(round.x[[i]])
	#image(dist.m[[i]])
	#plot(sort(dist.m[[i]]))
}