2010-11-10

多分岐を多次元で

R 木

こちらで２分岐木を描かせている
分岐木について考えてみる
枝の先で、２つに分かれれば２分岐木
numbif本に分かれれば、numbif-分岐木(たとえば、真菌)
２次元空間に広がるか(粘菌のように)か、３次元空間に広がるか(通常の木のように)：そういえば、もじほこりの乾燥胞子から寒天inタッパーで粘菌を育てていたこともあった。あのときの乾燥胞子は、今でも『生や』せるのだろうか…
多次元分岐のときに絵を描くとしたら、２次元に落とすしかない
２次元分岐木で、こちらでは、等長２次元回転行列と、短縮因子とでできた行列を用いている
２つに分岐するから、２つの行列を使っている
２次元での回転なので、２ｘ２行列を使っている
今、k次元で、numbif分岐させようと思ったら、kxk行列をnumbif個用意すればよい
また、描くときに、深度限定で、再帰関数を使うこともできる
以下のソースは、適当な次元空間に適当な分岐本数で適当な深度で描くためのもの
ただし、行列を作るのに乱数でいい加減な値を使っているので、「きれいに」描けていない。
きれいに描くためには、深度が進むにつれ、少し短くなる(0.8倍とか)ことと、分岐枝間の角度がそれなりであることが条件(その条件を入れるにはどうすればよいのかを放置したまま)
２次元、２分岐のときだけ、こちらで行列指定ができるように、特別処理をしている
掲載図は２次元、２分岐、深度８と３次元、５分岐、深度４

k<-3 # 空間次元
# 描くのはX,Yで指定する２次元
tobedrawnX<-1
tobedrawnY<-2
numbif<-5 # 枝分かれ本数
Niter<-4 # 枝分かれの深さ


library(MCMCpack)
lengthvector<-runif(numbif)*1/sqrt(k)
# 回転とは、長さkのベクトルを長さkのベクトルに移す行列のこと
BifMat<-array(0,c(k,k,numbif))
for(i in 1:numbif){
	#BifMat[,,i]<-rdirichlet(1,rep(1,k^2))*lengthvector[i]
	BifMat[,,i]<-rnorm(k^2)*lengthvector[i]
	# 次元k=2,枝分かれ本数numbif=2のときには、きれいに描けるように以下のようにしている
	if(k==2 && numbif==2){

		s<-pi/4    #右の枝の角度
		t<--pi/4*0.6  #左の枝の角度

		r<-0.8    #右の枝の長さ
		l<-0.7     #左の枝の長さ

		BifMat[,,1]<-r*matrix(c(cos(t),sin(t),-sin(t),cos(t)),2)
		BifMat[,,2]<-l*matrix(c(cos(s),sin(s),-sin(s),cos(s)),2)
	}
}



a<-b<-rep(0,k) # 第１、２点位置座標
b[1]<-1

# 描図範囲指定…まだ制御できていない
xlim<-ylim<-c(-5,5)
# 初めの２点をプロット
plot(a[tobedrawnX],a[tobedrawnY],xlim=xlim,ylim=ylim,cex=0.1)
par(new=TRUE)
plot(b[tobedrawnX],b[tobedrawnY],xlim=xlim,ylim=ylim,cex=0.1)
segments(a[tobedrawnX],a[tobedrawnY],b[tobedrawnX],b[tobedrawnY])

# 再帰的に点をプロットし、線でつないでいく
f<-function(a,b,itercounter,maxitercounter,numbif,k,BifMat,tobedrawnX,tobedrawnY,xlim,ylim){
	tmpa<-a
	tmpb<-b
	ret<-matrix(0,numbif,k)
	for(i in 1:numbif){
		ret[i,]<-tmpb+BifMat[,,i]%*%(tmpb-tmpa)
		#par(new=TRUE)
		#plot(tmpb[tobedrawnX],tmpb[tobedrawnY],xlim=xlim,ylim=ylim)
		par(new=TRUE)
		plot(ret[i,tobedrawnX],ret[i,tobedrawnY],xlim=xlim,ylim=ylim,cex=0.1)
		segments(tmpb[tobedrawnX],tmpb[tobedrawnY],ret[i,tobedrawnX],ret[i,tobedrawnY])
		if(itercounter+1<=maxitercounter){
			f(a=tmpb,b=ret[i,],itercounter=itercounter+1,maxitercounter=maxitercounter,numbif=numbif,k=k,BifMat=BifMat,tobedrawnX=tobedrawnX,tobedrawnY=tobedrawnY,xlim=xlim,ylim=ylim)
		}
		#print(ret[i,])
	}
	#ret

}

f(a=a,b=b,itercounter=1,maxitercounter=Niter,numbif=numbif,k=k,BifMat=BifMat,tobedrawnX=tobedrawnX,tobedrawnY=tobedrawnY,xlim=xlim,ylim=ylim)