せっかくなので、一番コンパクトな書き方をしてみよう。そして、少しづつ書き方を変えてみる。速度差は出るだろうか。
まず、そのまま書き下してみる。
計算環境はJulia Boxを使用
https://www.juliabox.com
し、計算機環境は
versioninfo()
Julia Version 0.6.0 Commit 9036443 (2017-06-19 13:05 UTC) Platform Info: OS: Linux (x86_64-pc-linux-gnu) CPU: Intel(R) Xeon(R) CPU E5-2673 v3 @ 2.40GHz WORD_SIZE: 64 BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell) LAPACK: libopenblas64_ LIBM: libopenlibm LLVM: libLLVM-3.9.1 (ORCJIT, haswell)
となっている。
二次元イジング模型のメトロポリス法によるモンテカルロシミュレーションのコードは、
function スイープ!(配置,Lx,Ly,w,mz)
for ix in 1:Lx
for iy in 1:Ly
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ifelse(ix==Lx,1,ix+1),iy]+配置[ifelse(ix==1,Lx,ix-1),iy]+配置[ix,ifelse(iy==Ly,1,iy+1)]+配置[ix,ifelse(iy==1,Ly,iy-1)])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
end
end
return 配置,mz
end
function モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
w = exp.(-[i for i in -8:8]/T) #ボルツマン重み e^[-ΔE/T]
mz = sum(配置)
Mz = 0
for i=1:熱化ステップ
配置,mz=スイープ!(配置,Lx,Ly,w,mz)
end
for i=1:最大ステップ
配置,mz=スイープ!(配置,Lx,Ly,w,mz)
Mz += abs(mz)
end
return Mz/(Lx*Ly*最大ステップ)
end
モンテカルロ法! (generic function with 1 method)
となり、このくらい短く書くことができる。
Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
106.624404 seconds (1.11 M allocations: 36.195 MiB, 0.01% gc time)
0.9920440608
次に、乱数を少し先に確保してみる。
function スイープ!(配置,Lx,Ly,w,mz)
r = rand(Lx,Ly) #乱数確保
for ix in 1:Lx
for iy in 1:Ly
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ifelse(ix==Lx,1,ix+1),iy]+配置[ifelse(ix==1,Lx,ix-1),iy]+配置[ix,ifelse(iy==Ly,1,iy+1)]+配置[ix,ifelse(iy==1,Ly,iy-1)])
配置[ix,iy] = ifelse(r[ix,iy] <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
end
end
return 配置,mz
end
スイープ! (generic function with 1 method)
Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
114.048818 seconds (3.03 M allocations: 74.626 GiB, 2.45% gc time)
0.988445197
少し遅くなった? 念のためもう一回実行する。
Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
115.313045 seconds (3.00 M allocations: 74.625 GiB, 2.53% gc time)
0.988445197
やはり遅くなっている。配列を確保したのが原因だろうか?
乱数に関しては元に戻して、次に、ifelse関数を三項演算子「a ? b : c 」に変えてみる
function スイープ!(配置,Lx,Ly,w,mz)
for ix in 1:Lx
for iy in 1:Ly
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ix==Lx?1:ix+1,iy]+配置[ix==1?Lx:ix-1,iy]+配置[ix,iy==Ly?1:iy+1]+配置[ix,iy==1?Ly:iy-1])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
end
end
return 配置,mz
end
スイープ! (generic function with 1 method)
Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
118.716406 seconds (1.01 M allocations: 31.153 MiB, 0.01% gc time)
0.9920440608
遅い気がする。念のためもう一度。
Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
121.250223 seconds (1.00 M allocations: 30.525 MiB, 0.01% gc time)
0.9920440608
うーん。これは遅い。つまり、現時点では
function スイープ!(配置,Lx,Ly,w,mz)
for ix in 1:Lx
for iy in 1:Ly
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ifelse(ix==Lx,1,ix+1),iy]+配置[ifelse(ix==1,Lx,ix-1),iy]+配置[ix,ifelse(iy==Ly,1,iy+1)]+配置[ix,ifelse(iy==1,Ly,iy-1)])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
end
end
return 配置,mz
end
スイープ! (generic function with 1 method)
Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
104.852433 seconds (1.01 M allocations: 31.152 MiB, 0.01% gc time)
0.9920440608
が一番速い。
次に、ifelse関数を余りを計算するmod1(ix-1,Lx)で置き換えてみよう。この関数は
ix = 1
println(mod1(ix-1,Lx))
ix = Lx
println(mod1(ix+1,Lx))
100 1
と周期境界条件を表現できる。これを使ってみる。
function スイープ!(配置,Lx,Ly,w,mz)
for ix in 1:Lx
for iy in 1:Ly
σi = 配置[ix,iy]
ΔE = 2σi*(配置[mod1(ix+1,Lx),iy]+配置[mod1(ix-1,Lx),iy]+配置[ix,mod1(iy+1,Ly)]+配置[ix,mod1(iy-1,Ly)])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
end
end
return 配置,mz
end
スイープ! (generic function with 1 method)
Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
628.313394 seconds (1.01 M allocations: 31.204 MiB, 0.00% gc time)
0.9920440608
とても遅くなった。ifelse関数が一番速い。
次に、コードが少し汚くなるが、ifelse関数も除去してみよう。
function スイープ!(配置,Lx,Ly,w,mz)
ix = 1
iy = 1
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,iy+1]+配置[ix,Ly])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
for iy in 2:Ly-1
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,iy+1]+配置[ix,iy-1])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
end
iy = Ly
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,1]+配置[ix,iy-1])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
for ix in 2:Lx-1
iy = 1
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ix+1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,Ly])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
for iy in 2:Ly-1
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ix+1,iy]+配置[iy-1,iy]+配置[ix,iy+1]+配置[ix,iy-1])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
end
iy = Ly
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ix+1,iy]+配置[ix-1,iy]+配置[ix,1]+配置[ix,iy-1])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
end
ix = Lx
iy = 1
σi = 配置[ix,iy]
ΔE = 2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,Ly])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
for iy in 2:Ly-1
σi = 配置[ix,iy]
ΔE = 2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,iy-1])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
end
iy = Ly
σi = 配置[ix,iy]
ΔE = 2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,1]+配置[ix,iy-1])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)
return 配置,mz
end
スイープ! (generic function with 1 method)
Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
101.889886 seconds (1.04 M allocations: 32.397 MiB, 0.01% gc time)
0.9990793904
すごい汚いコードになったが速くはなった。どうせならさらにいってみよう。mzの部分はもう少し高速化できる。
function スイープ!(配置,Lx,Ly,w,mz)
ix = 1
iy = 1
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,iy+1]+配置[ix,Ly])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += 配置[ix,iy]-σi
for iy in 2:Ly-1
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,iy+1]+配置[ix,iy-1])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += 配置[ix,iy]-σi
end
iy = Ly
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,1]+配置[ix,iy-1])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += 配置[ix,iy]-σi
for ix in 2:Lx-1
iy = 1
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ix+1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,Ly])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += 配置[ix,iy]-σi
for iy in 2:Ly-1
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ix+1,iy]+配置[iy-1,iy]+配置[ix,iy+1]+配置[ix,iy-1])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += 配置[ix,iy]-σi
end
iy = Ly
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ix+1,iy]+配置[ix-1,iy]+配置[ix,1]+配置[ix,iy-1])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += 配置[ix,iy]-σi
end
ix = Lx
iy = 1
σi = 配置[ix,iy]
ΔE = 2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,Ly])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += 配置[ix,iy]-σi
for iy in 2:Ly-1
σi = 配置[ix,iy]
ΔE = 2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,iy-1])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += 配置[ix,iy]-σi
end
iy = Ly
σi = 配置[ix,iy]
ΔE = 2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,1]+配置[ix,iy-1])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += 配置[ix,iy]-σi
return 配置,mz
end
スイープ! (generic function with 1 method)
Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
101.138176 seconds (1.04 M allocations: 32.277 MiB, 0.01% gc time)
0.9990793904
これは余り速くならなかった。しかし、あえてifelse関数を使う必要はないので、この変更は採用しよう。
結局、見やすさと速さを両立しているのは、
function スイープ!(配置,Lx,Ly,w,mz)
for ix in 1:Lx
for iy in 1:Ly
σi = 配置[ix,iy]
ΔE = 2σi*(配置[ifelse(ix==Lx,1,ix+1),iy]+配置[ifelse(ix==1,Lx,ix-1),iy]+配置[ix,ifelse(iy==Ly,1,iy+1)]+配置[ix,ifelse(iy==1,Ly,iy-1)])
配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)
mz += 配置[ix,iy]-σi
end
end
return 配置,mz
end
function モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
w = exp.(-[i for i in -8:8]/T) #ボルツマン重み e^[-ΔE/T]
mz = sum(配置)
Mz = 0
for i=1:熱化ステップ
配置,mz=スイープ!(配置,Lx,Ly,w,mz)
end
for i=1:最大ステップ
配置,mz=スイープ!(配置,Lx,Ly,w,mz)
Mz += abs(mz)
end
return Mz/(Lx*Ly*最大ステップ)
end
モンテカルロ法! (generic function with 1 method)
Lx = 100
Ly = 100
T = 1.0
最大ステップ = 1000000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
T = 1.0
@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
101.150472 seconds (1.02 M allocations: 31.624 MiB, 0.01% gc time)
0.9920440608
この形が良さそうだ。
using Plots
pyplot()
Plots.PyPlotBackend()
Lx = 100
Ly = 100
最大ステップ = 100000
熱化ステップ = 200
srand(123)
配置 = rand([-1,1],Lx,Ly)
nt = 20
time = [(nt+1-i)*0.2 for i in 1:nt]
Mztdep = zeros(Float64,nt)
i = 0
@time for T in time
i += 1
Mztdep[i] = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)
end
plot(time,Mztdep)