覆い尽くされるかどうかシミュレーション

• 昨日の記事で円盤に被覆される点の話を書いた
• 平面上の最密充填円配置の扱い、実際に乱点を発生させて、乱点の配置ごとに、最密充填円シートのランダム爆撃によるサバイバル個体数の分布がどうなるか(昨日の話は、無限に爆撃したら必ず…という設定だったが、有限回なら、乱点の配置のさせ方によって、生存個体数の期待値は変わるだろうし、最少生き残り個体数の分布も変わるだろう。そんなのをやるためのやっつけソース
• 固定した最密充填シートに対して、点の座標から最近接円盤の中心座標と距離を出す関数

nearestCenter <- function(X){
	scale.y <- X[2]/(sqrt(3))
	y.L <- floor(scale.y)
	y.U <- ceiling(scale.y)
	half.x <- X[1]/2
	half.x.2 <- (X[1]+1)/2
	tmpx <- round(half.x)*2
	tmpx2 <- round(half.x.2)*2-1
	ev <- FALSE
	if(y.L%%2 == 0){
		ev <- TRUE
	}
	if(ev){
		d1 <- (tmpx-X[1])^2+(y.L*sqrt(3)-X[2])^2
		d2 <- (tmpx2-X[1])^2+(y.U*sqrt(3)-X[2])^2
		if(d1>d2){
			ret <- c(tmpx2,y.U*sqrt(3))
			ret2 <- sqrt(d2)
		}else{
			ret <- c(tmpx,y.L*sqrt(3))
			ret2 <- sqrt(d1)
		}
	}else{
		d1 <- (tmpx-X[1])^2+(y.U*sqrt(3)-X[2])^2
		d2 <- (tmpx2-X[1])^2+(y.L*sqrt(3)-X[2])^2
		if(d1>d2){
			ret <- c(tmpx2,y.L*sqrt(3))
			ret2 <- sqrt(d2)
		}else{
			ret <- c(tmpx,y.U*sqrt(3))
			ret2 <- sqrt(d1)
		}
	}
	return(list(nearest = ret,dist=ret2))
}

• 上の関数がうまく行っているかの確認。乱点配置は、十分に広くしないと、被覆されやすくなる

n.pt <- 100000
L <- 1000 # 乱点の置ける範囲を指定するパラメタ
X <- matrix(runif(n.pt)*1000,ncol=2)

Z <- X
D <- rep(0,length(X[,1]))
for(i in 1:length(X[,1])){
	tmp <- nearestCenter(X[i,])
	Z[i,] <- tmp$nearest
	D[i] <- tmp$dist
}

plot(X)
for(i in 1:length(X[,1])){
	segments(X[i,1],X[i,2],Z[i,1],Z[i,2])
}
hist(D)

length(which(D>1))/length(D)

• 最密充填シートの爆撃は、シートのある点基準にして、その点を二次元についてランダムに動かすという操作と、その点についてランダムに回転するという操作の二つで、全パターンができるだろう(そうするとパラメタを３つ使うことになるのだが…自由度は２なので、どこか無駄がある…)

# 繰り返し回数。何度も乱点を発生させる
n.iter <- 100000
L <- 10000
survives <- rep(0,n.iter)
for(ii in 1:n.iter){
	n.pt <- 10 # 乱点数
	X <- matrix(runif(n.pt)*L,ncol=2)
	Z <- X
	D <- rep(0,length(X[,1]))
	for(i in 1:length(X[,1])){
		tmp <- nearestCenter(X[i,])
		Z[i,] <- tmp$nearest
		D[i] <- tmp$dist
	}
	# 生き残り人数を数える
	survives[ii] <- length(which(D>1))
}

hist(survives)

• どういう乱点配置がよいかは、幾度も配置してみて、それぞれが、多数回のシート爆撃に対して、どのようなサバイバル個体数分布をするかの集計が必要

# 乱点セットの発生回数
n.iter <- 1000
L <- 10
n.pt <- 20
# 乱点セットに対して行うランダム爆撃
n.rep <- 100
suv.list <- matrix(0,n.rep,n.pt+1)
X.hist <- list()
for(iii in 1:n.rep){

	X <- matrix(runif(n.pt*2)*L,ncol=2)
	X.hist[[iii]] <- X
survives <- rep(0,n.iter)
for(ii in 1:n.iter){
	x <- runif(2)*2
	theta <- runif(1)*2*pi
	R <- matrix(c(cos(theta),sin(theta),-sin(theta),cos(theta)),2,2)
	p <- t(t(X)-x)
	X.rot <- t(R %*% t(p))

	
	Z <- X.rot
	D <- rep(0,length(X[,1]))
	for(i in 1:length(X[,1])){
		tmp <- nearestCenter(X.rot[i,])
		Z[i,] <- tmp$nearest
		D[i] <- tmp$dist
	}
	
	survives[ii] <- length(which(D>1))
}
suv.list[iii,] <- tabulate(survives+1,n.pt+1)

hist(survives)
}
# 生存個体数の期待値
m <- apply(suv.list,1,function(x){sum(x*(0:(length(x)-1)))})
# 生存個体数の最少値
mmins <- apply(suv.list,1,min)
hist(m)
hist(mmins)

m.ord <- order(m,decreasing=TRUE)

plot(X.hist[[m.ord[1]]])
plot(X.hist[[m.ord[n.rep]]])