球の分裂と球面度数分布

分割

• 球の分裂に関するバナッハ=タルスキーパラドクスに端を発した球の分裂
• その新規の球の由来が元の球から均一にとられているのかどうかを評価するには、球面上で度数分布がとりたくなる
• そのために球面上の度数分布をとるときのタイルの配置の仕方とそのカウントの仕方、それを図示する方法について、書いた
• ここ数日の記事のこと
• それを併せる

library(rgl)
###########
# 2次元酔歩をさせよう
# 4種類の「一歩」があり得るので、それを作る
a <- 0.01
b <- 0.01
# 簡単のためにz軸を軸とした回転aと、y軸を軸とした回転bとを作る
# その回転行列
Ma <- matrix(c(cos(a),-sin(a),0,sin(a),cos(a),0,0,0,1),byrow=TRUE,3,3)
Mb <- matrix(c(cos(b),0,-sin(b),0,1,0,sin(b),0,cos(b)),byrow=TRUE,3,3)

# z軸の周りに正・逆２方向の回転を、y軸の周りにも２方向の回転をさせる
# a,-a,b,-bとする
# それぞのの回転行列をリストに納める

Ms <- list()
Ms[[1]] <- Ma
Ms[[2]] <- t(Ma)
Ms[[3]] <- Mb
Ms[[4]] <- t(Mb)

X0 <- runif(3)
X0 <- X0/sqrt(sum(X0^2))
Xall <- matrix(c(1,1,1,-1,-1,-1),byrow=TRUE,nrow=2)
s.list <- list()
cols <- rep(0,2)
n.iter <- 120
n.step <- 10000
for(ii in 1:n.iter){
	s <- sample(1:4,n.step,replace=TRUE)
	s[1] <- ii %% 4 +1
	tobeavoid <- c(2,1,4,3)
	s[2] <- sample((1:4)[-tobeavoid[s[1]]],1)
	s[3] <- sample((1:4)[-tobeavoid[s[2]]],1)
	
	# 回転軸を共有する正・逆回転が続くとそれは「不動」なので、そのような２歩を除く
	# この処理はしなくても、軌道的には問題ないが、「自由群」的には、ないことを持ち込む必要があるので、そのようにする
	s.2 <- c(s[1])
	for(i in 2:length(s)){
		tmp <- FALSE
		if(s[i-1] ==1){
			if(s[i] != 2){
				tmp <- TRUE
			}
		}else if(s[i-1] == 2){
			if(s[i] != 1){
				tmp <- TRUE
			}
		}else if(s[i-1] ==3){
			if(s[i] != 4){
				tmp <- TRUE
			}
		}else{
			if(s[i] != 3){
				tmp <- TRUE
			}
		}
		if(tmp) s.2 <- c(s.2,s[i])
	}
	length(s)-length(s.2)

	s <- s.2
	s.list[[ii]] <- s
	# 単位球表面にランダムに点を選び、そこから、酔歩させる
	X <- matrix(0,length(s)+1,3)

	X[1,] <- X0

	for(i in 2:length(X[,1])){
		X[i,] <- Ms[[s[i-1]]] %*% X[i-1,]
	}
	Xall <- rbind(Xall,X)
	cols <- c(cols,0,0,rep(s[1],length(s)+1-3),0)
}

plot3d(Xall,type="l",col=cols)


Sphere.Histogram.3.2 <- function(X,ts,k){
	ts.2 <- c(ts,ts[length(ts):1])
	cum.ts <- cumsum(ts)
	max.cum.ts <- cum.ts[length(cum.ts)]
	phis <- c(0,acos(1-cum.ts/max.cum.ts))
	diff.phis <- diff(phis)
	phis <- c(phis[-length(phis)],cumsum(c(phis[length(phis)],diff.phis[length(diff.phis):1])))
	#print(r)
	Hist <- list()
	
	thetas <- list()
	for(i in 1:k){
		thetas[[i]] <- seq(from = 0, to = 2*pi, length =ts[i]+1)
		Hist[[i]] <- rep(0,length(thetas[[i]]-1))
	}
	tmpks <- k:1
	for(i in 1:length(tmpks)){
		thetas[[k+i]] <- thetas[[k+1-i]]
		Hist[[k+i]] <- Hist[[k+1-i]]
	}
	print(thetas)
	Cmpx <- complex(real = X[,1],imaginary = X[,2])
	Md <- acos(X[,3])
	#Md <- Mod(Cmpx)
	Ag <- Arg(Cmpx)
	neg <- which(Ag < 0)
	Ag[neg] <- 2*pi +Ag[neg]
	for(i in 1:length(Hist)){
		s <- which(Md >= phis[i] & Md < phis[i+1])
		a <- Ag[s] %/% (thetas[[i]][2]) +1
		Hist[[i]] <- tabulate(a,nbins = ts.2[i]+1)
		#if(i<=k){
		#	Hist[[i]] <- tabulate(a,nbins = ts[i]+1)
		#}else{
		#	Hist[[i]] <- tabulate(a,nbins = ts[k+(ts[i]-k]+1)
		#}
		Hist[[i]] <- Hist[[i]][-length(Hist[[i]])]
	}
	return(list(Hist=Hist,X=X,Md=Md,Ag=Ag,t=t,k=k,phis=phis,thetas=thetas))
}
plot.SphereHist <- function(sp.out,ts,k,n=100000){
	n.pt <- n
	phis <- sp.out$phis
	thetas <- sp.out$thetas
	X <- matrix(rnorm(n.pt*3),ncol=3)
	X <- X/sqrt(apply(X^2,1,sum))
	Cmpx <- complex(real = X[,1],imaginary = X[,2])
	Md <- acos(X[,3])
	#Md <- Mod(Cmpx)
	Ag <- Arg(Cmpx)
	neg <- which(Ag < 0)
	Ag[neg] <- 2*pi +Ag[neg]
	cols <- rep(0,length(X[,1]))
	for(i in 1:length(sp.out$Hist)){
		s <- which(Md >= phis[i] & Md < phis[i+1])
		a <- Ag[s] %/% (thetas[[i]][2]) +1
		cols[s] <- sp.out$Hist[[i]][a]
	}
	library(rgl)
	plot3d(X,col = gray((max(cols)-cols)/(max(cols)-min(cols))),cex=0.1)
}


X <- Xall
t.len <- 5
ts <- c(1,6*(1:t.len))
k <- length(ts)
sp.out <- Sphere.Histogram.3.2(X,ts,k)

plot.SphereHist(sp.out,ts,k)

• 細胞分裂の映画(こちらやこちら)