# ただのメモ

library(sphere)

ThetaAngle<-function(cos,sin){
t<-acos(cos)
t[which(sin<0)]<-2*pi-t[which(sin<0)]
t
}

Simplex<-function(Ts,a,r){
A<-Ts%*%t(a)
R<-Ts%*%r
L<-sqrt(A^2+R^2)
Cs<-A/L
Ss<-R/L
Thetas<-ThetaAngle(Cs,Ss)
#Thetas[which(Ss<0)]<-2*pi -Thetas[which(Ss<0)]

TPointsX<-1/L * Cs
TPointsY<-1/L * Ss
xlim<-ylim<-c(-10,10)
plot(TPointsX,TPointsY,xlim=xlim,ylim=ylim)

selected<-c(which(A==max(A)))
current<-selected[1]
VsX<-c(1/max(A))
VsY<-c(0)
loop<-TRUE
while(loop){
NegativeStart<-FALSE

if(length(selected)==1 & Thetas[current]>pi){
NegativeStart<-TRUE
Mores<-c(which(Thetas>Thetas[current] & Thetas <(Thetas[current]+pi)), which(Thetas>=0 & Thetas <(Thetas[current]-pi)) )
}else{
NegativeStart<-FALSE
Mores<-which(Thetas>Thetas[current] & Thetas <(Thetas[current]+pi))
}
#if(Thetas[current]<pi){
#Mores<-which(Thetas>Thetas[current] & Thetas <(Thetas[current]+pi))
#}else{
#Mores<-c(which(Thetas>Thetas[current] & Thetas <(Thetas[current]+pi)), which(Thetas>=0 & Thetas <(Thetas[current]-pi)) )
#}

tmpL<-sqrt((1/L[Mores]*Cs[current]-1/L[current]*Cs[Mores])^2 +
(1/L[Mores]*Ss[current]-1/L[current]*Ss[Mores])^2 )

tmpSs<-(-1/L[current]*Cs[Mores]+1/L[Mores]*Cs[current]) / tmpL
tmpCs<-(1/L[current]*Ss[Mores]-1/L[Mores]*Ss[current]) / tmpL

tmpThetas<-ThetaAngle(tmpCs,tmpSs)
if(Thetas[current]+pi/2-2*pi>0){
tmpThetas[which(tmpThetas<Thetas[current]+pi/2-2*pi)]<-tmpThetas[which(tmpThetas<Thetas[current]+pi/2-2*pi)]+2*pi
}

#print(tmpThetas)
#print(NegativeStart)
#if(!NegativeStart){
#tmpThetas[which(tmpThetas<Thetas[current])]<-2*pi-tmpThetas[which(tmpThetas<Thetas[current])]
#tmpThetas[which(tmpThetas>pi]<-2*pi-tmpThetas[which(tmpThetas>pi]
#}

#if(Thetas[current]>pi){
#tmpThetas<-2*pi -tmpThetas
#}
#tmpThetas[which(tmpSs<0)]<-2*pi - tmpThetas[which(tmpSs<0)]

#tmpMoreThetas<-tmpThetas[which(tmpThetas>Thetas[current])]
if(length(tmpThetas)>0){
tmpid<-which(tmpThetas==min(tmpThetas))
tmpselected<-Mores[tmpid]
selected<-c(selected,tmpselected)
rx<-1/L[current]/cos(tmpThetas[tmpid]-Thetas[current])
VsX<-c(VsX,rx*cos(tmpThetas[tmpid]))
VsY<-c(VsY,rx*sin(tmpThetas[tmpid]))
current<-tmpselected

}else{
loop<-FALSE
}
if(current==selected[1]){
loop<-FALSE
}
#print(selected)
#plot(1/L[selected]*Cs[selected],1/L[selected]*Ss[selected])
for(i in 1:length(selected)){
abline(1/(L[selected[i]]*Ss[selected[i]]),-Cs[selected[i]]/Ss[selected[i]],col="red")
}
}
#print(Thetas[selected])

crossX<-crossY<-rep(0,length(selected)+1)
for(i in 1:length(selected)){
x1<-1/L[selected[i]]*Cs[selected[i]]
y1<-1/L[selected[i]]*Ss[selected[i]]
x2<-1/L[selected[i+1]]*Cs[selected[i+1]]
y2<-1/L[selected[i+1]]*Ss[selected[i+1]]
crossX[i]<-(y1*(x2^2+y2^2)-y2*(x1^2+y1^2))/(x2*y1-x1*y2)
crossY[i]<--(x1*(x2^2+y2^2)-x2*(x1^2+y1^2))/(x2*y1-x1*y2)
}
x1<-1/L[selected[length(selected)]]*Cs[selected[length(selected)]]
y1<-1/L[selected[length(selected)]]*Ss[selected[length(selected)]]
x2<-1/L[selected[1]]*Cs[selected[1]]
y2<-1/L[selected[1]]*Ss[selected[1]]
crossX[length(crossX)]<-(y1*(x2^2+y2^2)-y2*(x1^2+y1^2))/(x2*y1-x1*y2)
crossY[length(crossX)]<--(x1*(x2^2+y2^2)-x2*(x1^2+y1^2))/(x2*y1-x1*y2)

list(L=L,Cs=Cs,Ss=Ss,selected=selected,VsX=crossX,VsY=crossY)
}

#Simplex2
#最短辺を選ぶ
#ついで、現在の法線ベクトルから反時計周りにpiの法線ベクトルを持つ辺を対象とし
#交点と、現在の法線ベクトルとのなす余弦が最大になるような辺を次の辺として採用する
#最短辺に戻ったら終わり

Simplex2<-function(Ts,a,r){
A<-Ts%*%t(a)
R<-Ts%*%r
L<-sqrt(A^2+R^2)
Cs<-A/L
Ss<-R/L
Thetas<-ThetaAngle(Cs,Ss)
#Thetas[which(Ss<0)]<-2*pi -Thetas[which(Ss<0)]

TPointsX<-1/L * Cs
TPointsY<-1/L * Ss
xlim<-ylim<-c(-10,10)
plot(TPointsX,TPointsY,xlim=xlim,ylim=ylim)

#selected<-c(which(A==max(A)))
selected<-which(L==max(L))[1]
StartTheta<-Thetas[selected[1]]
Theta2<-(Thetas-StartTheta)%%(2*pi)
#Theta2[selected[1]]<-2*pi
current<-selected[1]
#VsX<-c(1/max(A))
#VsY<-c(0)
VsX<-c()
VsY<-c()
loop<-TRUE
#first<-TRUE
while(loop){
#if(!first){
#Theta2[selected[1]]<-2*pi
#}
#if(first){
#first<-FALSE
#}
currentTheta2<-Theta2[current]
#print(current)
plusPiTheta2<-currentTheta2+pi
EndTheta2<-min(2*pi,plusPiTheta2)
#print(EndTheta2)
if(EndTheta2==(2*pi)){ # 最初の法線からpi以上回ったら、最初の法線も考慮の対象に入れるため
Theta2[selected[1]]<-(currentTheta2+EndTheta2)/2
}
Mores<-which(Theta2>currentTheta2 & Theta2<EndTheta2)

tmpL<-sqrt((1/L[Mores]*Cs[current]-1/L[current]*Cs[Mores])^2 +
(1/L[Mores]*Ss[current]-1/L[current]*Ss[Mores])^2 )

tmpSs<-(-1/L[current]*Cs[Mores]+1/L[Mores]*Cs[current]) / tmpL
tmpCs<-(1/L[current]*Ss[Mores]-1/L[Mores]*Ss[current]) / tmpL

if(length(tmpCs)>0){
ttcc<-tmpCs * Cs[current] + tmpSs * Ss[current]
tmpid<-which(ttcc==max(ttcc))[1]
#tmpid<-which(tmpThetas==min(tmpThetas))
tmpselected<-Mores[tmpid]
selected<-c(selected,tmpselected)
current<-tmpselected

}else{
loop<-FALSE
}
if(current==selected[1]){
loop<-FALSE
}
#print(selected)
#plot(1/L[selected]*Cs[selected],1/L[selected]*Ss[selected])
#for(i in 1:length(selected)){
#abline(1/(L[selected[i]]*Ss[selected[i]]),-Cs[selected[i]]/Ss[selected[i]],col="red")
#}
}
#print(Thetas[selected])

crossX<-crossY<-rep(0,length(selected)-1)
for(i in 1:(length(selected)-1)){
x1<-1/L[selected[i]]*Cs[selected[i]]
y1<-1/L[selected[i]]*Ss[selected[i]]
#if(i<length(selected)){
x2<-1/L[selected[i+1]]*Cs[selected[i+1]]
y2<-1/L[selected[i+1]]*Ss[selected[i+1]]
#}else{
#x2<-1/L[selected[1]]*Cs[selected[1]]
#y2<-1/L[selected[1]]*Ss[selected[1]]
#}
crossX[i]<-(y1*(x2^2+y2^2)-y2*(x1^2+y1^2))/(x2*y1-x1*y2)
crossY[i]<--(x1*(x2^2+y2^2)-x2*(x1^2+y1^2))/(x2*y1-x1*y2)
}

list(L=L,Cs=Cs,Ss=Ss,selected=selected,VsX=crossX,VsY=crossY)
}

# 2次元凸包をある点から眺めた角度で考える
# 少しずつ角度を変えてそのときの、原点からの距離と
# 視点からの距離とを出し
# そのリストを返す

DistSimplex<-function(vs,p,n){# vs:(x,y)coordinates, p:new center, n: 各辺をn等分(角として)する
newvs<-t(t(vs)-p)
r<-sqrt(apply(newvs^2,1,sum))
#outx<-outy<-rep(0,n*length(vs[,1]))
outx<-outy<-matrix(0,length(vs[,1]),n)
counter<-0
deltatheta<-rep(0,length(vs[,1]))
for(i in 1:length(vs[,1])){
#for(i in 1:18){
st<-i
end<-i+1
if(i==length(vs[,1])){
end<-1
}
if(r[st]!=0 & r[end]!=0){
cosst<-newvs[st,1]/r[st]
sinst<-newvs[st,2]/r[st]
cosend<-newvs[end,1]/r[end]
sinend<-newvs[end,2]/r[end]

tst<-ThetaAngle(cosst,sinst)
tend<-ThetaAngle(cosend,sinend)

difft<-tend-tst
if(difft>=0 & difft <=pi){

}else if(difft > pi & difft<=2*pi){
difft<-(-1)*(2*pi - difft)
}else if(difft<0 & difft >=-pi){

}else if(difft< -pi & difft >= -2*pi){
difft<-2*pi+difft
}
deltatheta[i]<-difft
#if(difft<0){
#difft<-difft+2*pi
#}
#print(difft)
ts<-tst+difft*(0:(n-1))/n + difft/(2*n)
newt1<-tst-ts
newt2<-tend-ts
newrs<-r[st]*r[end]/(r[st]*sin(newt1)-r[end]*sin(newt2)) * sin(difft)
#if(i!=0){
#par(new=TRUE)
#}
#plot(newrs*cos(ts),newrs*sin(ts))
#outx[(counter+1):(n*i)]<-newrs*cos(ts)
#outy[(counter+1):(n*i)]<-newrs*sin(ts)
outx[i,]<-newrs*cos(ts)
outy[i,]<-newrs*sin(ts)
}
counter<-n*i
#plot(outx,outy)
}
return(list(x=outx+p[1],y=outy+p[2],p=p,n=n,rp=sqrt(outx^2+outy^2),deltatheta=deltatheta))
}

df<-5

Nt<-8

# 法線ベクトルの集合
Ts<-RandomSphere(n=Nt,df=df)
Ts<-rbind(Ts,-Ts)

# 真実のモデルの方向の単位ベクトル
a<-RandomSphere(n=1,df=df)

# 輪切りにする面の数
Ns<-3
# その単位ベクトル
rs<-RandomSphere(n=Ns,df=df)
ip<-rs %*% t(a)
rs<-t(t(rs)-t(a)%*%t(ip))

rs<-rs / sqrt(apply(rs^2,1,sum))

plot(as.data.frame(rs))
# 統計量の大きさとモデルの強さとの比として計算する数(Nr*2+1)
Nr<-8

Ratios<-10^((-Nr):Nr)

# 統計量とモデルの強さの比あたり、計算する、角度の数
Na<-12
Angles<-seq(from=0,to=1,by=1/(Na-1))*2*pi

for(i in 1:length(rs[,1])){

simplexOut<-Simplex2(Ts,a,rs[i,])

#simplexOut<-Simplex(Ts,a,r)

#simplexOut
xlim<-ylim<-c(min(1/simplexOut$L*simplexOut$Cs,1/simplexOut$L*simplexOut$Ss),
max(1/simplexOut$L*simplexOut$Cs,1/simplexOut$L*simplexOut$Ss))

xlim<-ylim-c(min(simplexOut$VsX,simplexOut$VsY),max(simplexOut$VsX,simplexOut$VsY))

xlim<-ylim<-c(-4,4)
plot(1/simplexOut$L*simplexOut$Cs,1/simplexOut$L*simplexOut$Ss,xlim=xlim,ylim=ylim,pch=20)
for(i in 1:length(simplexOut$selected)){ abline(1/(simplexOut$L[simplexOut$selected[i]]*simplexOut$Ss[simplexOut$selected[i]]),-simplexOut$Cs[simplexOut$selected[i]]/simplexOut$Ss[simplexOut$selected[i]],col="red") } par(new=TRUE) plot(simplexOut$VsX,simplexOut$VsY,xlim=xlim,ylim=ylim,pch=19,col="blue") } plot(simplexOut$VsX,simplexOut$VsY,pch=19,col="blue") for(i in 1:Ns){ #単位凸包を作る simplexOut<-Simplex2(Ts,a,rs[i,]) #print(UnitSimplex) xlim<-ylim<-c(-10,10) plot(1/simplexOut$L*simplexOut$Cs,1/simplexOut$L*simplexOut$Ss,xlim=xlim,ylim=ylim,pch=20) for(i in 1:length(simplexOut$selected)){
abline(1/(simplexOut$L[simplexOut$selected[i]]*simplexOut$Ss[simplexOut$selected[i]]),-simplexOut$Cs[simplexOut$selected[i]]/simplexOut$Ss[simplexOut$selected[i]],col="red")
}

par(new=TRUE)
plot(simplexOut$VsX,simplexOut$VsY,xlim=xlim,ylim=ylim,pch=19,col="blue")

#相対的な位置を変えたとき(Ratiosの値)の関数として表したい
#アウトプットは、２交点の座標なのだが・・・虚数も含めて-b \pm sqrt(D) /2a のDとかを使って表す

}

################
# A/K ratiosを出す
# 等比級数にすることで、計算量を減らす

geomSeries<-function(stval,endval,n){
exp(seq(from=log(stval),to=log(endval),length.out=n))
}
maxp<-10^(-1)
minp<-10^(-15)

minchival<-sqrt(qchisq(maxp,df=df,lower.tail=FALSE))
maxchival<-sqrt(qchisq(minp,df=df,lower.tail=FALSE))
n<-10

geomSer<-geomSeries(minchival,maxchival,n)

logchisqval<-seq(from=log(minchisqval),to=log(maxchisqval),length.out=10)

alphachival<-exp(logchisqval)
betachival<-c(0,alphachival)

# betachival/alphachival の値は、alphachivalとbetachivalが同じ等比数列なので、比の-(要素数)乗・・・(要素数)乗になる

ratio<-alphachival[1]/alphachival[2]

ratios<-ratio^((-length(alphachival)+1):(length(alphachival)-1))

outprob<-matrix(0,length(ratios),length(betachival))

for(x in 1:Ns){
simplexOut<-Simplex2(Ts,a,rs[x,])
# 凸包の形は決まった
vs<-cbind(simplexOut$VsX,simplexOut$VsY)

#area<-rep(0,length(ratios))
for(i in 1:length(ratios)){
p<-c(ratios[i],0)

n<-100

# 視点からの等角点の計算
dsout<-DistSimplex(vs,p,n)
xlim<-ylim<-c(min(dsout$x,dsout$y,p),max(dsout$x,dsout$y,p))
plot(c(t(dsout$x)),c(t(dsout$y)),type="l",xlim=xlim,ylim=ylim)
par(new=TRUE)
plot(p[1],p[2],col="red",xlim=xlim,ylim=ylim)

#area[i]<-sum((pchisq(dsout$rp^2,df=df))*dsout$deltatheta/(2*pi*dsout$n)) for(j in 1:length(outprob[1,])){ outprob[i,j]<-outprob[i,j]+sum((pchisq((dsout$rp*betachival[j])^2,df=df))*dsout$deltatheta/(2*pi*dsout$n))

}
}
matplot(outprob,type="l")

}
outprob<-outprob/Ns

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