■

• 数値列がある。それを入れ替えたい
• すべての要素が入れ替えられていて、一つとしてもとの順番になっていないようにする
  • 使える関数があるかと思ったが見つからないので、作る
  • バグがあったので少し修正しました(20110126)

shuffle<-function(x){# xは数値ベクトル
	X<-length(x) # 要素数
	s<-sample(1:X) # １から要素数の数列
	ret<-1:X #入れ替えた順番を納めるベクトル
	evod<-X%%2 # 要素数の偶奇
	L<-(X-evod)/2 # ペアワイズに入れ替えるので、ペアの数をとる
	for(i in 1:L){# 入れ替えた順列の先頭から順にペアを作って交換する
		pre<-s[i*2-1]
		post<-s[i*2]
		ret[pre]<-post
		ret[post]<-pre
	}
	if(evod==1){ 
# 要素数が奇数のときはペアになれないので、その要素だけ、最後にどれでもよいので入れ替える
		post<-s[X]
		ds<-1:(X-1)
		t<-sample(ds[-which(ret==X)],1)
		pre<-ret[s[t]]
		ret[s[t]]<-post
		ret[post]<-pre
	}
	x[ret]
}
# 確かめる
x<-1:9
# たくさんshuffle()関数を適用する
Niter<-10000
cp<-rep(0,Niter)
for(i in 1:Niter){
	sx<-shuffle(x)
# 前後の値の差の最小値を試行ごとに格納する
	cp[i]<-min(abs(x-sx))
# 途中でxとsxとで同じ値になることがあったら、知らせる(こうならない)
	if(cumprod(x-sx)[length(x)]==0)print(sx)

}
# shuffle()関数がうまくいっていることをshuffle()前後値が０にはならないことで確認する
min(abs(cp))

• この作業は任意次元の正規直交基底を回転させるため
• さらにそれは、曲線にあって接線方向の運動方向を決めるため
• 以下のソースはほんのメモ

####################
shuffle<-function(x){
	X<-length(x)
	s<-sample(1:X)
	ret<-1:X
	evod<-X%%2
	L<-(X-evod)/2
	for(i in 1:L){
		pre<-s[i*2-1]
		post<-s[i*2]
		ret[pre]<-post
		ret[post]<-pre
	}
	if(evod==1){
		t<-sample(s[-X],1)
		pre<-ret[t]
		post<-ret[X]
		ret[t]<-post
		ret[X]<-pre
	}
	x[ret]
}

x<-1:10

Niter<-10000
cp<-rep(0,Niter)
for(i in 1:Niter){
	sx<-shuffle(x)

	cp[i]<-min(abs(x-sx))
	if(cumprod(x-sx)[length(x)]==0)print(sx)

}
min(abs(cp))
Niter<-10000
dt<-0.01
d<-2
xs<-matrix(0,Niter,d)

xs[1,]<-runif(d)

n<-3
ctrs<-matrix(runif(d*n),ncol=d)

M<-list()

for(i in 1:n){
	ok<-TRUE
	news<-rep(0,d)
	while(ok){
		s<-1:d
		news<-rep(0,d)
		already<-c()
		for(j in 1:(d-1)){
			already2<-c(already,j)
			print(already2)
			print(s[-already2])
			news[j]<-sample(s[-already2],1)
			already<-c(already,news[j])
		}
		if(s[-already]!=d){
			news[d]<-s[-already]
			ok<-FALSE
		}
	}
	print(news)
	cumprod(news-s)
	pm<-1
	M[[i]]<-matrix(0,d,d)
	for(j in 1:d){
		M[[i]][j,news[j]]<-pm
		pm<-pm*(-1)
	}
}

M


for(i in 2:Niter){
	tmp<-rep(0,d)
	for(j in 1:n){
		tmpr<-sqrt(sum(xs[i-1,]-ctrs[j,])^2)

		tmpv<-M[[j]]%*%(xs[i-1,]-ctrs[j,])/sqrt(sum(tmpv^2))
		print(sum(tmpv*(xs[i-1,]-ctrs[j,])))
		tmp<-tmp+tmpv/tmpr*dt
	}
	xs[i,]<-xs[i-1,]+dt*tmp
}

xs2<-rbind(xs,ctrs)
plot(as.data.frame(xs2))