多様体推定

• トーラスの一部(少し曲がった円筒)の３次最小全域木を作ってみた
• 処理はすごく重い：重くて使えない
• 多様体推定がどんな具合なのかはだいたいわかったので、いい方法がないか、探してみることにする
  • たとえばこちら
• Isomapとは→こちら
• RでIsomap→veganパッケージのisomap()こちら

SimplexVolume<-function(x,Factorial=TRUE){
	n<-length(x[,1])
	#d<-t(x[2:n,])-x[1,]
	d<-apply(x,2,FUN="diff")
	if(Factorial){
		ret<-log(abs(det(d))) - lfactorial(n-1)
	}else{
		ret<-log(abs(det(d)))
	}
	return(ret)
}


SimplexVolume2<-function(x,Factorial=TRUE){
	n<-length(x[,1])
	tmpX<-t(x[2:n,])-x[1,]
	if(n==2){
		ret<-log(sqrt(sum(tmpX^2)))
	}else{
		tmpX<-SubDfSVD(t(tmpX))
		#svd.out<-svd(t(tmpX))
		#print(svd.out)
		#print(tmpX)
		#print(svd.out[[3]])
		#tmpX<-t(tmpX)%*%svd.out[[3]]
		#print(tmpX)

		ret<-SimplexVolume(rbind(rep(0,length(tmpX[1,])),tmpX),Factorial=Factorial)
	}
	return(ret)
}

SubDfSVD<-function(RR){
# 特異値分解
	svd.out<-svd(t(RR))
# 特異値分解の結果と、点の元の座標を用いて
# k次元座標にする
	RR%*%svd.out[[2]]
}



library(sets)
library(rgl)
# No. nodes
Nv<-10
# Dimensions
df<-9
# Coordinates
#X<-matrix(runif(Nv*df),ncol=df)
X<-matrix(rnorm(df*Nv),ncol=df)
nrm<-sqrt(apply(X^2,1,sum))
X<-X/nrm
#t<-runif(Nv)*2*pi
#X<-cbind(cos(t),sin(t))
X[,1]<-abs(X[,1])
X[,2]<-abs(X[,2])
EX<-matrix(runif(Nv*df),ncol=df)
EX<-CategoryVector(Nv)
EX<-X+matrix(rnorm(ncol(X)*nrow(X)),ncol=ncol(X))*1
df<-ncol(X)


##########
library(GPArotation)
Nv<-10
df<-10
X<-matrix(0,Nv,df)
kkk<-4
for(i in 1:kkk){
	X[,i]<-runif(Nv)
}
#X[,1]<-runif(Nv)
#X[,2]<-runif(Nv)
plot3d(X[,1:3])

RM<-Random.Start(df)
X<-X%*%RM
##########
R1<-10
R2<-3

Nv<-200

X<-matrix(0,Nv,3)
for(i in 1:Nv){
	t1<-runif(1)*2*pi
	t2<-runif(1)*2*pi
	tmpx<-R1*cos(t1)
	tmpy<-R1*sin(t1)
	tmpz<-R2*cos(t2)
	tmpx<-tmpx+R2*sin(t2)*cos(t1)
	tmpy<-tmpy+R2*sin(t2)*sin(t1)
	X[i,]<-c(tmpx,tmpy,tmpz)
}
plot3d(X)


################
library(sets)
library(rgl)
library(gtools)
library(GPArotation)

Nv<-30
t<-seq(from=0,to=1,length=Nv)*2*pi*10

x<-sin(t)+rnorm(length(t))*0.1
y<-cos(t)+rnorm(length(t))*0.1
z<-rep(0,length(x))


plot(x,y)

X<-cbind(x,y,z)
df<-3

RM<-Random.Start(df)
X<-X%*%RM
Nv<-length(x)

##########
R1<-10
R2<-1
df<-3
Nv<-200

X<-matrix(0,Nv,3)
for(i in 1:Nv){
	t1<-runif(1)*2*pi/12
	t2<-runif(1)*2*pi
	tmpx<-R1*cos(t1)
	tmpy<-R1*sin(t1)
	tmpz<-R2*cos(t2)
	tmpx<-tmpx+R2*sin(t2)*cos(t1)
	tmpy<-tmpy+R2*sin(t2)*sin(t1)
	X[i,]<-c(tmpx,tmpy,tmpz)
}
plot3d(X)
RM<-Random.Start(df)
X<-X%*%RM
plot3d(X)


##############



Ngs<-list()

Factorial<-FALSE

GList<-list()
GVolume<-list()
GVolumeSum<-list()
k<-1

GList[[k]]<-list()
GVolume[[k]]<-list()
#GVolumeSum[[k]]<-0
for(i in 1:Nv){
	GList[[k]][[i]]<-as.set(c(i))
	GVolume[[k]][[i]]<-0
	#GVolumeSum[[k]]<-GVolumeSum[[k]]+GVolume[[k]][[i]]
}

ks<-1:df
for(kk in 1:length(ks)){
	k<-ks[kk]
	Ngs[[k]]<-length(GList[[k]])
	#print(GList)
	#print(GList[[k]])

	currentGList<-GList[[k]]
	currentNgs<-length(currentGList)
	#Em<-Ws<-matrix(0,Ngs[[k]],Ngs[[k]])
	Em<-Ws<-matrix(Inf,currentNgs,currentNgs)
	#for(i in 1:(Ngs[[k]]-1)){
	for(i in 1:(currentNgs-1)){
		#for(j in (i+1):Ngs[[k]]){
		for(j in (i+1):currentNgs){
			#dif1<-GList[[k]][[i]]-GList[[k]][[j]]
			#dif2<-GList[[k]][[j]]-GList[[k]][[i]]
			dif1<-currentGList[[i]]-currentGList[[j]]
			dif2<-currentGList[[j]]-currentGList[[i]]
			if(length(dif1)==1 & length(dif2)==1){
				Em[i,j]<-1
				#vs<-set_union(GList[[k]][[i]],GList[[k]][[j]])
				vs<-set_union(currentGList[[i]],currentGList[[j]])
				Ws[i,j]<-Ws[j,i]<-SimplexVolume2(X[unlist(vs),],Factorial=Factorial)
				
			}
		}
	}
	#diag(Ws)<-max(Ws)+1
	#print(Em)
	#print(Ws)

	cnt<-1
	GList[[k+1]]<-list()
	GVolume[[k+1]]<-list()
	#GVolumeSum[[k+1]]<-0
	Candidates<-rep(1,Ngs[[k]])
	NonCandidates<-rep(0,Ngs[[k]])

	for(i in 1:(Ngs[[k]]-1)){
		#print(currentGList)
		currentNgsX<-currentNgs
		currentNgs<-length(currentGList)
		NewCandidates<-rep(0,currentNgs)
		NewNonCandidates<-rep(1,currentNgs)
		NewCandidates[1:currentNgsX]<-Candidates
		NewNonCandidates[1:currentNgsX]<-NonCandidates
		Candidates<-NewCandidates
		NonCandidates<-NewNonCandidates
		#print(Candidates)
		#print(NonCandidates)
		#print(length(which(Candidates==1)))
		#print(length(which(NonCandidates==1)))
		#print("---")
		#print(currentNgsX)
		#print(currentNgs)
		#Em<-Ws<-matrix(0,Ngs[[k]],Ngs[[k]])
		if(i!=1){
			WsX<-Ws
			Em<-Ws<-matrix(Inf,currentNgs,currentNgs)
			Ws[1:nrow(WsX),1:ncol(WsX)]<-WsX
			if(currentNgs>currentNgsX){
				#for(i in 1:(Ngs[[k]]-1)){
				#for(ii in 1:(currentNgs-1)){
				for(ii in which(Candidates==1)){
					#for(j in (i+1):Ngs[[k]]){
					for(jj in max(currentNgsX+1,ii+1):currentNgs){
						#dif1<-GList[[k]][[i]]-GList[[k]][[j]]
						#dif2<-GList[[k]][[j]]-GList[[k]][[i]]
						dif1<-currentGList[[ii]]-currentGList[[jj]]
						dif2<-currentGList[[jj]]-currentGList[[ii]]
						if(length(dif1)==1 & length(dif2)==1){
							Em[ii,jj]<-1
							#vs<-set_union(GList[[k]][[i]],GList[[k]][[j]])
							vs<-set_union(currentGList[[ii]],currentGList[[jj]])
							Ws[ii,jj]<-Ws[jj,ii]<-SimplexVolume2(X[unlist(vs),],Factorial=Factorial)
							
						}
					}
				}
			}
			

		}else{
			Em<-Ws<-matrix(Inf,currentNgs,currentNgs)
			#for(i in 1:(Ngs[[k]]-1)){
			for(ii in 1:(currentNgs-1)){
				#for(j in (i+1):Ngs[[k]]){
				for(jj in (ii+1):currentNgs){
					#dif1<-GList[[k]][[i]]-GList[[k]][[j]]
					#dif2<-GList[[k]][[j]]-GList[[k]][[i]]
					dif1<-currentGList[[ii]]-currentGList[[jj]]
					dif2<-currentGList[[jj]]-currentGList[[ii]]
					if(length(dif1)==1 & length(dif2)==1){
						Em[ii,jj]<-1
						#vs<-set_union(GList[[k]][[i]],GList[[k]][[j]])
						vs<-set_union(currentGList[[ii]],currentGList[[jj]])
						Ws[ii,jj]<-Ws[jj,ii]<-SimplexVolume2(X[unlist(vs),],Factorial=Factorial)
						
					}
				}
			}
		}
		#diag(Ws)<-max(Ws)+1
		#Ws[which(Ws==0)]<-max(Ws)+1
		#print(Em)
		#print(Ws)

		if(i==1){
			tmpVolume<-min(Ws[upper.tri(Ws)])
			GVolume[[k]][[i]]<-exp(tmpVolume)
			#GVolumeSum[[k]]<-GVolumeSum[[k]]+GVolume[[k]][[i]]
			selected<-which(Ws==tmpVolume,arr.ind=TRUE)
			if(length(selected)>2)selected<-selected[1,]
			Candidates[selected]<-0
			NonCandidates[selected]<-1
			#print(Candidates)
			#print(NonCandidates)
		}else{
			tmpWs<-Ws[which(Candidates==1),which(NonCandidates==1)]
			#print(tmpWs)
			tmpVolume<-min(tmpWs)
			tmpselected<-which(tmpWs==tmpVolume,arr.ind=TRUE)
			#print(tmpselected)
			if(length(tmpselected)>2)tmpselected<-tmpselected[1,]
			GVolume[[k]][[i]]<-exp(tmpVolume)
			#GVolumeSum[[k]]<-GVolumeSum[[k]]+GVolume[[k]][[i]]
			if(i!=(Ngs[[k]]-1)){
				selected<-c(which(Candidates==1)[tmpselected[1]],which(NonCandidates==1)[tmpselected[2]])
			}else{
				selected<-c(which(Candidates==1),which(NonCandidates==1)[tmpselected[1]])
			}
			
			#print(tmpselected)
			#print(selected)
			Candidates[selected[1]]<-0
			NonCandidates[selected[1]]<-1
			#print(Candidates)
			#print(NonCandidates)
		}
		#print(GList[[k]][[selected[1]]])
		#print(GList[[k]][[selected[2]]])
		#GList[[k+1]][[cnt]]<-set_union(GList[[k]][[selected[1]]],GList[[k]][[selected[2]]])
		GList[[k+1]][[cnt]]<-set_union(currentGList[[selected[1]]],currentGList[[selected[2]]])
		#Ws[selected]<-max(Ws[upper.tri(Ws)])+1
		# エッジを足すと、k次元の完全グラフがk-1個増える
		#tmpset<-set_intersection(GList[[k]][[selected[1]]],GList[[k]][[selected[2]]])
		tmpset<-set_intersection(currentGList[[selected[1]]],currentGList[[selected[2]]])
		currentGListLen<-length(currentGList)
		tmpcnt<-1
		for(j in tmpset){
			currentGList[[currentGListLen+tmpcnt]]<-GList[[k+1]][[cnt]]-j
			tmpcnt<-tmpcnt+1
		}
		cnt<-cnt+1

	}
	GVolumeSum[[k]]<-sum(unlist(GVolume[[k]]))

	#print(GList[[k+1]])
}

GList

plot(X[,1],X[,2])

for(i in length(GList):2){
	for(j in 1:length(GList[[i]])){
		scmb<-set_combn(GList[[i]][[j]],2)
		#print(scmb)
		for(j2 in scmb){
			#print(j2)
			v<-c()
			for(j3 in j2){
				v<-c(v,j3)
			}
			segments(X[v[1],1],X[v[1],2],X[v[2],1],X[v[2],2],col=i)
		}
	}
}

plot3d(X[,1],X[,2],X[,3])
for(i in length(GList):2){
	for(j in 1:length(GList[[i]])){
		scmb<-set_combn(GList[[i]][[j]],2)
		#print(scmb)
		for(j2 in scmb){
			#print(j2)
			v<-c()
			for(j3 in j2){
				v<-c(v,j3)
			}
			segments3d(X[v,],col=i)
		}
	}
}

GVolumeSum

incr<-unlist(GVolumeSum)[2:length(GVolumeSum)]/unlist(GVolumeSum)[1:(length(GVolumeSum)-1)]

par(mfcol=c(2,2))
plot(incr)
plot(unlist(GVolumeSum))
plot(GVolumeSum[[1]]^(1:(length(GVolumeSum)-1)))
plot(unlist(GVolume))
par(mfcol=c(1,1))