周期データを考えるときのいろいろ

• 昨日の続き
• 多次元球座標の関数と多次元トーラスの関数を別々に書いた
• 統一して書くと以下のようになるらしい
• k次元トーラス、k次元球の場合には、「角」の変数がk-1個必要なのは、トーラスも球も同じで、座標の増やし方の部分で

ret<-k*ret+C[i]*incr

• というように、「これまでのret」に増分を加えていく形になっている。
• この「これまでのret」をどの程度重みづけするかという観点から言うと、球は重みづけなし、トーラスは重みづけありという違い(らしい)
• ただし、球と「球のように見えるトーラス」とでは、本質的に異なるので、曲線はつながらずに「ワープ」するようだ→「r<-0.8 # これを1にすると、連続性がなくなる(ようだ)」

• 少し書き換え(20110126)

# トーラスの前に多次元球を描こう
# k次元球の座標はオイラーの多次元角座標を用いて、k-1個の角パラメタと
# １個の半径パラメタとの合わせてk個のパラメタで表せる
# 単位k次元球の場合は、半径が１なので、
# 角座標(k-1個の角)ｄ表される
# このk次元球面の直交座標系の座標を出す関数は以下の通り
# vはk-1個の角度のベクトル

sphereCoords<-function(v){
	# すべての角度のコサインとサインをまず出す
	C<-cos(v)
	S<-sin(v)
	# xk=sin(vk)
	# xk-1=cos(vk)*sin(vk-1)
	# xk-2=cos(vk)*cos(vk-1)*sin(vk-2)
	# ...
	# x2=cos(vk)*cos(vk-1)*cos(vk-2)*...*sin(v1)
	# x1=cos(vk)*cos(vk-1)*cos(vk-2)*...*cos(v1)
	# なので以下の通り
	ret<-c(1,S)

	for(i in 1:length(v)){
		for(j in i:length(v)){
			ret[i]<-ret[i]*C[j]
		}
	}
	return(ret)
}

# これを用いて、(k次元)トーラスの座標を角座標で作る
# まず、ドーナツの主要な輪について、回転角を与える
# 次に、主要な輪の点を中心として、その位置ベクトルを含む平面であって
# 輪を乗せている面と垂直な平面を考え
# そこにドーナツの実の部分の輪を描く
# 第２の角度はこの実の部分の輪における角度である
# ただし、トーラスの場合には、輪の階層ごとに異なる半径があるのに対し
# 多次元球の場合は、角座標を上げて行っても半径は同じであることに注意
# 普通のドーナツの場合はここまで
# 次元を上げるときは、作業を繰り返す

# トーラスk次

# パラメタは以下の通り
# いわゆるトーラス
# 周期の重層性の度数
k<-2 # kを増やしてもよいようにパラメタ化してあるが、ドーナツはk=2

# トーラスでは、主要な輪について処理した後、
# その位置を中心として次の輪の分の座標を足しこむが、
# そうしなければいけないというわけでもないので、
# 輪の強さをだんだん小さくするための係数を incrRatioとする
incrRatio<-1
# トーラス上の直線は、ドーナツの主要な輪をある周波数で回りつつ、ドーナツの実の部分の小さな輪もある周波数で回る
# その周波数を与えるためのパラメタがPs
Ps<-c(7,11)

# ドットを打つ時点を作る
n<-1000
t<-seq(from=0,to=1,length.out=n)*2*pi



MultDimTorus3<-function(v,k){
	C<-cos(v)
	S<-sin(v)
	ret<-rep(0,length(v)+1)
	preret<-ret
	for(i in 1:length(v)){
		if(i==1){
			ret[1]<-C[1]
			ret[2]<-S[1]
		}else{
			# 次元を上げながら座標を作り上げていくときに
			# 前段階に増えて分の座標を確認する
			# 新たな輪の分は、次のルールで足しこむ
			# 前段階までに使われていない軸には新な輪のコサインを与える
			# 残りの足しこみ分はサインだが、それは、前段階までに使われた軸に
			# 分配する必要がある
			# その分配が「前段階に増えた分の座標」に応じて行われる
			incr<-ret-preret 
			preret<-ret
			incr<-incr/sqrt(sum(incr^2))
			ret<-k*ret+C[i]*incr
			ret[i+1]<-S[i]
			
		}
	}
	return(ret)
}

# パラメタを振って描いてみよう
# k=2で固定する
k<-2
incrRatio<-4
# 主要な輪は３回まわって元に戻る
# ドーナツの実は５回まわって元に戻る
Ps<-c(3,5) 

n<-10000

t<-seq(from=0,to=10,length.out=n)*2*pi


Xs<-matrix(0,length(t),k+1)
for(i in 1:length(t)){
	Xs[i,]<-MultDimTorus3(t[i]*Ps,incrRatio)
}
xlim<-ylim<-zlim<-c(min(Xs),max(Xs))

plot3d(Xs[,1],Xs[,2],Xs[,3],cex=0.1,col=gray((1:length(t))/length(t)),xlim=xlim,ylim=ylim,zlim=zlim)

plot(as.data.frame(Xs))
matplot(Xs,type="l")

# 普通にトーラス
incrRatio<-4
# 主要な輪は３回まわって元に戻る
# ２つの輪の周期が無理数だと、元に戻らないので、ぐるぐる回る
# 無理数としてpiとexpをとる
Ps<-c(pi,exp(1)) 
t<-seq(from=0,to=10,length.out=n)*2*pi

Xs<-matrix(0,length(t),k+1)
for(i in 1:length(t)){
	Xs[i,]<-MultDimTorus3(t[i]*Ps,incrRatio)
}
xlim<-ylim<-zlim<-c(min(Xs),max(Xs))

plot3d(Xs[,1],Xs[,2],Xs[,3],cex=0.1,col=gray((1:length(t))/length(t)),xlim=xlim,ylim=ylim,zlim=zlim)

plot(as.data.frame(Xs))
matplot(Xs,type="l")

# incrRatioを1にすると、ドーナツの中心がなくなる
k<-2
incrRatio<-1
Ps<-c(pi,exp(1)) 
t<-seq(from=0,to=10,length.out=n)*2*pi

Xs<-matrix(0,length(t),k+1)
for(i in 1:length(t)){
	Xs[i,]<-MultDimTorus3(t[i]*Ps,incrRatio)
}
xlim<-ylim<-zlim<-c(min(Xs),max(Xs))

plot3d(Xs[,1],Xs[,2],Xs[,3],cex=0.1,col=gray((1:length(t))/length(t)),xlim=xlim,ylim=ylim,zlim=zlim)

plot(as.data.frame(Xs))
matplot(Xs,type="l")

# incrRatioを0にすると、球の表面のみを通過する
k<-2
incrRatio<-0
Ps<-c(pi,exp(1)) 
t<-seq(from=0,to=10,length.out=n)*2*pi

Xs<-matrix(0,length(t),k+1)
for(i in 1:length(t)){
	Xs[i,]<-MultDimTorus3(t[i]*Ps,incrRatio)
}
xlim<-ylim<-zlim<-c(min(Xs),max(Xs))

plot3d(Xs[,1],Xs[,2],Xs[,3],cex=0.1,col=gray((1:length(t))/length(t)),xlim=xlim,ylim=ylim,zlim=zlim)

plot(as.data.frame(Xs))
matplot(Xs,type="l")