JuliaとRとの速さの比較をしてみる - ryamadaのコンピュータ・数学メモ

Dirichlet問題と言うのがある
周辺の値が固定されているときに内部で調和する関数はどんなものか、という問題
これをシミュレーションで解いてみることにする(正方形という簡単な領域での単純な例)
以下は、こちらにあるJuliaのソース(Juliaのインストールなどはこちら)と、それを翻訳したRのソース
どういうシミュレーションか、というと：
- 正方形領域の周囲に固定値を与える
- 正方形の内部の点については、
  - そこから酔歩をさせて、初めて行き着いた周辺の値が何かを何度も記録し、その平均を持って、そこの値とする
- どういう意味かというと、固定値を持つ周辺から、正方形内部の点に対して平等に影響が及ぶとする。その平等な影響の及び方が酔歩であるとする
このシミュレーションがRのベクトル処理とかでやったりするのに向かないのは、「酔歩がいつ辺縁に到着するかがわからないから、それを目途にできないため」
さて、十分に平均がとれるまで回数を回すとするとN=5000とかになるが、それだとjuliaは終わるがRで終わらないので、N=100でやってみるとすると
juliaが0.4秒くらい
Rが18秒強

#
# load within julia using include("dirichlet.jl")
#
function rytest(N)
nx, ny = 25, 25
f = zeros(nx,ny)
f[1,:], f[nx,:], f[:,1], f[:,ny] = 2, 2, -2, -2
C = zeros(nx,ny)
#N = 1
#
for x = 1:nx
 for y = 1:ny
  for k = 1:N
   xx = x
   yy = y
   while ((xx > 1) && (xx < nx) && (yy > 1) && (yy < ny))
    b = rand(2)
    if (b[1] < .5)
     if (b[2] < .5)
      xx -= 1
     else
      xx += 1
     end
    elseif (b[2] < .5)
     yy -= 1
    else
     yy += 1
    end
   end # while
   C[x,y] += f[xx,yy]
  end # for k
 end # for y
end # for x

draw(PNG("dirichlet.png", 10cm, 10cm), spy(C/N))
return 0
end


@time rytest(100)

rytest <- function(N){
nx <- ny <- 25
f <- matrix(0,nx,ny)
f[1,] <- f[nx,] <- 2
f[,1] <- f[,ny] <- -2
z <- matrix(0,nx,ny)
#N <- 100

for(i in 1:nx){
	for(j in 1:ny){
		for(k in 1:N){
			xx <- i
			yy <- j
			while((xx>1) & (xx<nx) & (yy>1) & (yy<ny)){
				r <- runif(2)
				if(r[1]>0.5){
					if(r[2]>0.5){
						xx <- xx + 1
					}else{
						xx <- xx -1
					}
				}else{
					if(r[2]>0.5){
						yy <- yy + 1
					}else{
						yy <- yy -1
					}
				}
			}
			z[i,j] <- f[xx,yy]
		}
	}
}

image(z/N)
return(0)
}

system.time(rytest(100))

１次元GFFのjuliaシミュレーションとそのR翻訳版
- 左右端が0になっていて、内部がランダム

#
# load within julia using include("gff.jl")
#
Pkg.add("Gadfly") # produces "INFO: Nothing to be done." if already installed.
Pkg.add("Cairo") # idem.
#Pkg.add("PyPlot") # idem.
using Compose, Gadfly, Cairo # for 2D graphics
#using PyPlot # for 3D graphics (Gadfly does not provide 3D graphics)
#
# dimension 1
#
L = 500
Z = randn(L-1)
lambda = 2*sin(pi*[1:L-1]/(2*L)).^2
evec = sqrt(2)*sin(pi*kron([1:L-1],[1:L-1]')/L)
GFF1 = zeros(L-1)
for x = 1:L-1
    for y = 1:L-1
        for n = 1:L-1
            GFF1[x] += (lambda[n])^(-1/2)*evec[n,x]*evec[n,y]*Z[y]
        end
    end
end
# graphics ith Gadfly
p = plot(x = [1:L-1], y = GFF1,
         Geom.line,
         Guide.xlabel(""),
         Guide.ylabel(""),
         #Guide.title("GFF"),
         Theme(default_color=color("black")))
draw(PNG("gff1.png", 10cm, 10cm), p)

L <- 500
Z <- rnorm(L-1)

xy <- expand.grid(1:(L-1),1:(L-1))
xy2 <- matrix(xy[,1]*xy[,2],L-1,L-1)

lambda <- 2*sin(pi*(1:(L-1))/(2*L)) ^2
evec <- sqrt(2)*sin(pi*xy[,1]*xy[,2]/L)
evec2 <- matrix(evec,L-1,L-1)

GFF1 <- rep(0,L-1)

for(i in 1:(L-1)){
	for(j in 1:(L-1)){
		for(k in 1:(L-1)){
			GFF1[i] <- GFF1[i] + lambda[k] ^ (-1/2)*evec2[k,i]*evec2[k,j]*Z[j]
		}
	}
}

plot(1:(L-1),GFF1,type="l")



library(rgl)
plot3d(xy[,1],xy[,2],z)