ループ盛りだくさんの処理を比較する

• こちらにself-avoiding pathのシミュレーションがある
• Rcppを使って比較してみる

#include <RcppArmadilloExtensions/sample.h>
// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp ;

// [[Rcpp::export]]

using namespace Rcpp;

IntegerMatrix selfAvoidingPathCpp2(int k,int past,int t,int maxt) {
	IntegerMatrix x(t,k);
	for(int i=0;i<k;++i){
		for(int j=0;j<t;++j){
			x(j,i) = 0;
		}
		
	}
	RNGScope scope;
	

	IntegerVector kv = seq_len(k);
	IntegerVector s(2);
	s[0] = -1;
	s[1] = 1;
	IntegerVector r = RcppArmadillo::sample(kv, maxt, TRUE);
	IntegerVector rs = RcppArmadillo::sample(s, maxt, TRUE);
	int cnt = 0;
	for(int i=0;i<maxt;++i){
		if(cnt==t-1){
			break;
		}
		IntegerVector here(k);
		for(int j=0;j<k;++j){
			here[j] = x(cnt,j);
		}
		IntegerVector toward(k);
		for(int j=0;j<k;++j){
			toward[j] = here[j];
		}
		toward[r[i]-1] = toward[r[i]-1]+rs[i];
		if(past == 0){
			++cnt;
			for(int j=0;j<k;++j){
				x(cnt,j) = toward[j];
			}
		}else{
			int tmppast = 0;
			if(past == -1){
				
			}else{
				int tmp2 =cnt-past;
				if(tmppast<tmp2){
					tmppast = tmp2;
				}
			}
			bool b = TRUE;
			for(int j=tmppast;j<cnt+1;++j){
				bool b2 = TRUE;
				for(int jj=0;jj<k;++jj){
					if(x(j,jj) != toward[jj]){
						b2 = FALSE;
						break;
					}
				}
				if(b2){
					b = FALSE;
					break;
				}
			}
			if(b){
				++cnt;
				for(int j=0;j<k;++j){
					x(cnt,j) = toward[j];
				}
			}else{
				++cnt;
				for(int j=0;j<k;++j){
					x(cnt,j) = here[j];
				}
			}
		}
		
	}


	return x;
}

library(Rcpp)
library(RcppArmadillo)
library(rgl)

sourceCpp("selfAvoidingPathCpp2.cpp")

• R版

# k:次元,past:考慮する過去(-1は無限に),t:歩数,max.t:最大トライ数
selfAvoidingWalkR <- function(k,past=-1,t=100,max.t=1000){
	X <- matrix(0,t,k)
	R <- sample(1:k,max.t,replace=TRUE)
	R.sign <- sample(c(-1,1),max.t,replace=TRUE)
	cnt <- 1
	for(i in 1:(max.t-1)){
		if(cnt==t){
			break
		}
		here <- X[cnt,]
		toward <- X[cnt,]
		toward[R[i]] <- toward[R[i]] + R.sign[i]
		if(past == 0){
			cnt <- cnt+1
			X <- X[cnt,] <- toward
		}else{
			if(past == -1){
				tmp.past <- 1
			}else{
				tmp.past <- cnt-past+1
				tmp.past <- max(1,tmp.past)
			}
			if(cnt==tmp.past){
				tmp.X <- matrix(X[cnt,],nrow=1)
			}else{
				tmp.X <- X[tmp.past:cnt,]
			}
			tmp.d <- apply(abs(t(tmp.X)-toward),2,sum)
			if(prod(tmp.d)!=0){
				cnt <- cnt+1
				X[cnt,] <- toward
			}else{
				cnt <- cnt+1
				X[cnt,] <- here
			}
		}
	}
	X
}

• • cppの方では、既通過点があるかないかの処理のところで、既通過点を順に調べて、あったらループを終了している(本当は近い過去から調べる方が効率的だが、今のソースは遠い過去から調べているので、さらに高速になるはず)のに対して、対象となるすべての過去の点を調べているところは、処理数の違いを生んでいるので、「対等な比較」ではないのだが、おそらくRの方では、既通過点を逐一調べると遅くなるだろうと予想するので、ここではそれは比較しないこととする)

> k<-2
> max.t<-10000
> past <- -1
> t <- 1000
> benchmark(
+ selfAvoidingPathCpp2(k,past,t,max.t),
+ selfAvoidingWalkR(k,past=past,t=t,max.t=max.t),
+ columns=c("test", "replications","elapsed", "relative"),
+ replications=100,
+ order="relative")
                                                     test replications
1                 selfAvoidingPathCpp2(k, past, t, max.t)          100
2 selfAvoidingWalkR(k, past = past, t = t, max.t = max.t)          100
  elapsed relative
1    0.06        1
2  100.68     1678