2010-05-03

多次元で垂線を下ろす２

N次元空間を考える
ある点Aがある
いくつかの点の集合Xsがある
- Xsが線形独立でN個のとき
- Xsの頂点を結んだ面へのAからの垂線の距離を求めるのが、"MinDistanceX()"
- N次元空間にN個の点があって、その頂点を結ぶ凸包をAと原点を結ぶ線が通過するかどうかの判定をするのが"JudgeInside()"
いくつかの点があって、それらのすべてのペアを結び一番外側をとったとき凸包を作る
- Aがこの凸包の表面か外にあるとき、凸包への最短距離を求めるのが"MinDistanceRecursiveX2()"(これは"MinDistanceX()"を使って、垂線を下ろしたときに内部に落ちればそれを距離とし、外部に落ちれば、下ろした先の「外周」に対して垂線を下ろして、一番近いところを探す、再帰関数)
内部の点からの距離は０
外部の点からの距離は有限
凸包上の点をランダムに発生させて、そこへのAからの距離が、求めた距離より大きいことも確かめる。

#R05.R
# 多次元空間において、凸包内にあるかどうかの判定をする
# 点の数が次元以上ならば、次元個ずつの点の組み合わせについて内部判定をする

JudgeInside<-function(A,Xs,bothside=TRUE){

 B<-NULL
 df<-length(A) # 次元
 N<-length(Xs[,1]) # 点の数
 inside<-FALSE
 coefs<-NULL
 if(N>=df){
  library(gtools)
  cmb<-combinations(N,df,repeats=FALSE)

  loop<-TRUE
  counter<-1
  vs<-NULL
  firstB<-TRUE
  while(counter<=length(cmb[,1])&loop){
#選んだ点の組み合わせを座標としてAの位置を表す
   tmp<-solve(t(Xs[cmb[counter,],]),A)
   if(firstB){
    B<-t(Xs[cmb[counter,],]) %*% tmp
   }
   if(length(tmp[tmp<0])==0 || (bothside & length(tmp[tmp>0])==0) ){# すべての要素が0以上の条件
    inside<-TRUE # 錐の内側
    vs<-Xs[cmb[counter,],]
    coefs<-tmp
    loop<-FALSE
   }
   counter<-counter+1
  }
  if(inside){
   B<-t(vs) %*% coefs
  }

 }else{
  tmp<-solve(t(Xs[,1:N]),A[1:N])

  if(sum((A-B)^2)<10^(-16)& (length(tmp[tmp<0])==0|| (bothside & length(tmp[tmp>0])==0))){
   inside=TRUE
   B<-t(Xs) %*% tmp
   coefs<-tmp
  }

 }
 list(inside=inside,A=A,Xs=Xs,coefs=coefs,B=B)
}

####examples
df<-6 # 次元
N<-4 # 点の数
 library(MCMCpack)
Xs<-matrix(rdirichlet(N,rep(1,df)),nrow=N)
#Xs<-matrix(c(1,0,0,1),nrow=2,byrow=TRUE)

if(N>=df){
#点の数が次元以上ならば、内側になるかどうかを試すべくAを作る
 A<-rnorm(df)
}else{
#点の数が次元未満なら、あえて、線形結合になるようなAを作ってやる(すべてが正なら内側)
 A<-apply(Xs *rnorm(N),2,sum)
}


JudgeInside(A,Xs)

################
#k次元空間中のk個のベクトルの作る凸包への垂線とその長さを求める

MinDistance2<-function(A,Xs,bothside=TRUE){
 df<-length(A)
 N<-length(Xs[,1])
 Xs1<-t(t(Xs)-A)
 C<-apply(Xs1,2,sum)/N
 Xs2<-t(t(Xs1)-C)
 Xs3<-Xs2 %*% t(Xs1)
 Xs4<-rbind(Xs3[1:(N-1),],rep(1,N))
 coefs<-solve(Xs4,c(rep(0,(N-1)),1))
 Apr<-t(Xs1) %*% coefs
 L<-nrm(Apr)
 Apr<-Apr+A
 inside<-FALSE
 if(length(coefs[coefs<0])==0 || (bothside & length(coefs[coefs>0])==0) ){# すべての要素が0以上の条件
  inside<-TRUE
 }
 list(A=A,Xs=Xs,distance=L,Apr=Apr,coefs=coefs,inside=inside)

}

#次元と点数が等しい場合のみ
MinDistance<-function(A,Xs,bothside=TRUE){
 N<-length(Xs[,1])
 C<-apply(Xs,2,sum)/N # Xsの重心
 Xs2<-t(t(Xs)-C)

# Xsのベクトルの頂点を通る面の法線ベクトルを求める
#Tansの最後の要素値を１に固定し、N-1連立方程式
 Tans<-c(solve(Xs2[1:(N-1),1:(N-1)],-Xs2[1:(N-1),N]),1)
 Tans<-Tans/sqrt(sum(Tans^2))
#Aの足Aprを求める。足の長さも
 tmpXmat<-as.matrix(cbind(Xs,rep(1,N)))

 tmpApr<-solve(as.matrix(cbind(t(tmpXmat),c(Tans,0))),c(A,1))
 k<-tmpApr[N+1]
 Apr<-t(Xs)%*% tmpApr[1:N]
 coefs<-tmpApr[1:N]
 inside<-FALSE
 if(length(coefs[coefs<0])==0 || (bothside & length(coefs[coefs>0])==0) ){# すべての要素が0以上の条件
  inside<-TRUE
 }
 L<-abs(k)
 list(A=A,Xs=Xs,distance=L,Apr=Apr,coefs=coefs,inside=inside,Tvector=Tans)
}

###examples
#k次元k線形独立ベクトルセット
#１個の外点
library(MCMCpack)
df<-N<-4 # 次元
#N<-df-2 # 点の数
#各行が点を表すベクトル
Xs<-matrix(runif(N*df),nrow=N)
#Xs<-matrix(rdirichlet(N,rep(1,df)),nrow=N)
A<-rnorm(df)

#df<-N<-2
#Xs<-matrix(runif(4),nrow=2)

MDout<- MinDistance(A,Xs)
MDout2<- MinDistance2(A,Xs)

MDout
MDout2

#Tansは面上の任意の点を結んだベクトルと直交する
tmptmprs<-rdirichlet(2,rep(1,length(A)))
tmptmprs<-matrix(runif(2*length(A)),nrow=2)
tmptmpsum<-apply(tmptmprs,1,sum)
tmptmprs<-tmptmprs/tmptmpsum

B1<-t(Xs)%*%tmptmprs[1,]
B2<-t(Xs)%*%tmptmprs[2,]
sum((B1-B2)*MDout$Tvector)


#求めた垂線の足が、AからXs面への最短距離であることを確認
df<-N<-4 # 次元
N<-df-2 # 点の数
#各行が点を表すベクトル
Xs<-matrix(runif(N*df),nrow=N)
#Xs<-matrix(rdirichlet(N,rep(1,df)),nrow=N)
A<-rnorm(df)
L<-MDout2(A,Xs,bothside=FALSE)$distance
testN<-1000
Ls<-rep(0,testN)
tmprs<-matrix(rnorm(testN*N)*100,nrow=testN)
tmpsum<-apply(tmprs,1,sum)
tmprs<-tmprs/tmpsum
for(i in 1:testN){
 tmpA<-t(Xs) %*% tmprs[i,]
 Ls[i]<-sqrt(sum((A-tmpA)^2))
}

plot(sort(Ls),ylim=c(min(L,Ls),max(L,Ls)),type="l")
abline(h=L,col="red")
 min(Ls)-L

#######################################
#単体への最短距離を再帰的に求める
#k次元にk+1個の点があり
#線形独立であること
#それらが単体(凸包)をなす
################
MinDistRecursive<-function(A,Xs,bothside=TRUE){
 library(gtools)
 df<-length(A)
 N<-length(Xs[,1])
 if(N==1){
  L<-sqrt(sum((Xs-A)^2))
  return(list(A=A,Xs=Xs,distance=L,Apr=Xs,coefs=c(1),inside=TRUE))
 }
 MDout<-MinDistance2(A,Xs,bothside)
 print(MDout$distance)
 if(MDout$inside){
  return(list(A=MDout$A,Xs=MDout$Xs,distance=MDout$distance,Apr=MDout$Apr,coefs=MDout$coefs,inside=MDout$inside))
 }else{
  select<-NULL

  for(i in 1:N){
   tmpout<-MinDistRecursive(A,matrix(Xs[which((1:N)!=i),],nrow=N-1),bothside)
   print(tmpout$distance)
   if(i==1){
    select<-tmpout
   }else{
    if(select$distance>tmpout$distance){
     select<-tmpout
    }
   }
  }
  return(list(A=select$A,Xs=select$Xs,distance=select$distance,Apr=select$Apr,coefs=select$coefs,inside=select$inside))
 }
}


#############
df<-N<-6 # 次元
N<-3 # 点の数
 library(MCMCpack)
Xs<-matrix(rdirichlet(N,rep(1,df)),nrow=N)
#Xs<-matrix(c(1,0,0,1),nrow=2,byrow=TRUE)
A<-rnorm(df)*100
#finalanswer2<-MinDistRecursive2(A,Xs)

finalanswer<-MinDistRecursive(A,Xs)

finalanswer

L<-finalanswer$distance


finalanswer$distance
#finalanswer2$distance



testN<-1000
Ls<-rep(0,testN)
#tmprs<-matrix(runif(testN*length(Xs[,1]))*100,nrow=testN)
tmprs<-rdirichlet(testN,rep(1,length(Xs[,1])))
#tmpsum<-apply(tmprs,1,sum)
#tmprs<-tmprs/tmpsum
for(i in 1:testN){
 tmpA<-t(Xs) %*% tmprs[i,]
 Ls[i]<-sqrt(sum((A-tmpA)^2))
}
plot(sort(Ls),ylim=c(min(L,Ls),max(L,Ls)))
abline(h=L,col="red")
min(Ls)-L