たとえば

• ソースは、非常に重いのだが…
```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)
}
}
}

```