たとえば

• ソースは、非常に重いのだが…

SimplexVolume<-function(x,Factorial=FALSE){
	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=FALSE){
	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)
# No. nodes
Nv<-20
# Dimensions
df<-3
# 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])
Ngs<-list()

GList<-list()

k<-1

GList[[k]]<-list()
for(i in 1:Nv){
	GList[[k]][[i]]<-as.set(c(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(0,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),])
				
			}
		}
	}
	diag(Ws)<-max(Ws)+1
	#print(Em)
	#print(Ws)

	cnt<-1
	GList[[k+1]]<-list()

	Candidates<-rep(1,Ngs[[k]])
	NonCandidates<-rep(0,Ngs[[k]])

	for(i in 1:(Ngs[[k]]-1)){
		#print(currentGList)
		currentNgsX<-currentNgs
		currentNgs<-length(currentGList)
		#print(currentNgsX)
		#print(currentNgs)
		#Em<-Ws<-matrix(0,Ngs[[k]],Ngs[[k]])
		if(i!=1){
			WsX<-Ws
			Em<-Ws<-matrix(0,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(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),])
							
						}
					}
				}
			}
			

		}else{
			Em<-Ws<-matrix(0,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),])
						
					}
				}
			}
		}
		#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),])
					
				}
			}
		}
		diag(Ws)<-max(Ws)+1
		#print(Em)
		#print(Ws)

		if(i==1){
			selected<-which(Ws==min(Ws[upper.tri(Ws)]),arr.ind=TRUE)
			Candidates[selected]<-0
			NonCandidates[selected]<-1
			#print(Candidates)
			#print(NonCandidates)
		}else{
			tmpWs<-Ws[which(Candidates==1),which(NonCandidates==1)]
			#print(tmpWs)
			tmpselected<-which(tmpWs==min(tmpWs),arr.ind=TRUE)
			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

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

}


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)
		}
	}
}