2011-02-25

小区画の総合成果

R 時系列解析周期トーラス干渉

生物の特長の一つは、「小部屋」を作ること
- 個体は臓器で構成される
- 臓器は組織で構成される
- 組織は細胞で構成される
- 細胞は小器官で構成される
- ...
- 個体は集団に属する
- 集団は大集団に属する
- 集団は他種集団とともに生態系に属する
どこの階層でもよいけれど、構成員が同じ「仕組み」を持っているとする
「仕組み」がうまく回るためには、「周期的」にしておくのは、良い方法だ(そうでないと、「静的な定常状態」に落ち着きやすいけれど、「(単純な)周期状態」は(比較的)作り出しやすいから(反応系をサイクルにして一方通行にする、とか)
構成要素の「仕組み」がそれぞれ回るときには、「仕組み」が持つ「基本周期」があるだろう
ただし、個々の構成要素がその「仕組み」をまわすときには、それぞれの構成要素の置かれた事情によって、「基本周期」にぴったり合わせることは難しく、ずれが生じることだろう
そのずれを一生懸命、「中央集権的」に制御することは、いかにも生物「らしくない」
放任しつつ、結果とし目的を果たさせたい
「構成要素」が「仕組み」をまわすか否かのON/OFFをすることで、並列で動く「仕組み」の個数を制御するのはどうだろう
「構成要素」が独立するように「小部屋に分ける」として、その分け方を階層的にすることも、実際に生物ではよくやられていることだ
さて。そんなことを実行するソース
- 掲載図は
  - 全体が５部屋に分かれ、５部屋の内部がそれぞれ10部屋に分かれたもの
  - 一番小さな単位の小部屋では、ある周期とその1/300くらいでランダムにずれを生じさせ
  - 各小部屋の周期の位相はまったくランダムにしたもので
  - 一番上の階層の濃度変化をプロットしたもの
- これを作るのに、「任意進数」を使用(こちらにソース)
こちらとこちらとも関連

# 木を根から生やそう
counterplus<-function(x,t){
	x[1]<-x[1]+1
	for(i in 1:(length(x))){
		if(x[i]==t[i]){
			x[i]<-0
			if(i<length(x)){
				x[i+1]<-x[i+1]+1
			}
			
		}else{
			return(x)
		}
	}
	return(x)
}
counterplus2<-function(x,t){
	x[1]<-x[1]+1
	for(i in 1:(length(x))){
		if(x[i]==(t[i]+1)){
			x[i]<-1
			if(i<length(x)){
				x[i+1]<-x[i+1]+1
			}
			
		}else{
			return(x)
		}
	}
	return(x)
}
# これが、階層的「仕切り」を定めるパラメタ
# 全体を５部屋に仕切って、次に10部屋に仕切って
# その小部屋に1種類の反応が回っている
nchild<-c(5,10,1)
nstep<-length(nchild)

counter<-rep(0,nstep)
loop<-TRUE
cnt<-1
ns<-NULL
ns[[cnt]]<-counter
while(loop){
	counter<-counterplus2(counter,nchild)
	print(counter)
	cnt<-cnt+1
	ns[[cnt]]<-c(counter)

	if(sum(counter)==sum(nchild)){
		loop<-FALSE
	}
}
cnt
ns
layer<-rep(0,length(ns))
# 周期を定める
w<-runif(1)*3
#d<-rnorm(nchild[nstep])
#w<-runif(length(ns))
# 部屋ごとの周期ブレ
d<-runif(length(ns))*0.1
# 上位階層に束ねるときに、小部屋の寄与率はランダムにする
r<-runif(length(ns))
# 階層的部屋(大部屋から小部屋まで)を行列形式のグラフ表現に
OyakoM<-matrix(0,length(ns),length(ns))
for(i in 1:(length(ns)-1)){
	for(j in (i+1):length(ns)){
		tmp<-ns[[i]]-ns[[j]]
		nonzero<-which(tmp!=0)
		if(length(nonzero)==1){
			if(length(which(tmp[(nonzero+1):nstep]==0))==nstep-nonzero){
				if(ns[[i]][nonzero]*ns[[j]][nonzero]==0){
					if(tmp[nonzero]>0){
						OyakoM[i,j]<--r[i]
						OyakoM[j,i]<-r[i]
					}else{
						OyakoM[i,j]<-r[j]
						OyakoM[j,i]<--r[j]
					}
				}
				
			}
		}
	}
}
# 時系列計算
L<-10000
cycle<-100
t<-seq(from=0,to=1,length.out=L)*2*pi*cycle

X<-NULL
for(i in 1:length(ns)){
	X[[i]]<-rep(0,length(t))
}
for(i in length(ns):1){
	#print(i)
	if(ns[[i]][nstep]==0){
	
	}else{
		tmpphase<-runif(1)*2*pi
		#print(tmpphase)
		#X[[i]]<-cos((w+d[ns[[i]][nstep]])*(t+tmpphase))
		X[[i]]<-cos((w+d[i])*(t+tmpphase))
	#	print(X[[i]])
	}
	for(j in 1:length(OyakoM[,i])){
		if(OyakoM[j,i]>0){
			#print(j)
			X[[j]]<-X[[j]]+OyakoM[j,i]*X[[i]]
		}
		
	}
}
xmat<-NULL
for(i in 1:length(ns)){
	xmat<-rbind(xmat,X[[i]])
}
for(i in 1:length(ns)){
	tmp<-length(which(ns[[i]]==0))
	layer[i]<-2^tmp
	xmat[i,]<-xmat[i,]/layer[i]
}

xmat<-t(xmat)
#matplot(xmat,type="l")


plot(xmat[,1],type="l")