多分岐を多次元で２

• 昨日の続き
• 分岐木を描くときに、伸びてきた枝から、ある回転をk次元空間で起こさないといけない
  • 回転とは、長さkのベクトルを長さkのベクトルに移す行列のこと
  • 正規直交基底行列が回転を表す
  • 正規直交基底行列は、k要素が作るk(k-1)/2ペアに関する回転を合わせたものとして表せることを使う

k<-7 # 次元
# 単位行列を作る
a<-matrix(0,k,k)
diag(a)<-1
a
# 適当にk(k-1)/2ペアに関する回転角を定める
matt<-matrix(runif(k*k),k,k)
# k(k-1)/2回転させる
for(i in 1:(k-1)){
	for(j in (i+1):k){
		t<-matt[i,j]
		rotm<-matrix(0,k,k)
		diag(rotm)<-1
		rotm[i,i]<-rotm[j,j]<-cos(t)
		rotm[i,j]<-(-sin(t))
		rotm[j,i]<-sin(t)
		#print(rotm)
		a<-rotm%*%a
		#print(a)
	}
}
a
# 単位ベクトルであることの確認
apply(a^2,1,sum)
apply(a^2,2,sum)
# 直交の確認
for(i in 1:(k-1)){
	for(j in (i+1):k){
		print(sum(a[i,]*a[j,]))
	}
}

• これを使って、昨日のむちゃくちゃな絵を、それなりに見せられるようにする

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


library(MCMCpack)
lengthvector<-runif(numbif)*0.2+0.8

# 回転とは、長さkのベクトルを長さkのベクトルに移す行列のこと
# 正規直交基底行列が回転を表す
# 正規直交基底行列は、k要素が作るk(k-1)/2ペアに関する回転を合わせたものとして表せることを使う

BifMat<-array(0,c(k,k,numbif))

for(x in 1:numbif){
	A<-matrix(0,k,k)
	diag(A)<-1
	#a
	matt<-matrix(runif(k*k),k,k)
	for(i in 1:(k-1)){
		for(j in (i+1):k){
			t<-matt[i,j]
			rotm<-matrix(0,k,k)
			diag(rotm)<-1
			rotm[i,i]<-rotm[j,j]<-cos(t)
			rotm[i,j]<-(-sin(t))
			rotm[j,i]<-sin(t)
			print(rotm)
			A<-rotm%*%A
			print(A)
		}
	}
	
	BifMat[,,x]<-lengthvector[x]*A
}

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(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,])
	}
}

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