versioninfo()
Julia Version 1.1.1 Commit 55e36cc308 (2019-05-16 04:10 UTC) Platform Info: OS: Windows (x86_64-w64-mingw32) CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz WORD_SIZE: 64 LIBM: libopenlibm LLVM: libLLVM-6.0.1 (ORCJIT, haswell) Environment: JULIA_CMDSTAN_HOME = C:\CmdStan JULIA_NUM_THREADS = 4 JULIA_PKGDIR = C:\JuliaPkg
# 2D Ising Metropolis (Gen Kuroki, 2017-01-09)
const β_crit = log(1+sqrt(2))/2
rand_ising2d(n) = rand([-1, 1], n, n)
function ising2d_sum_of_adjacent_spins(s, m, n, i, j)
s[ifelse(i+1 ≤ m, i+1, 1), j] + s[ifelse(i-1 ≥ 1, i-1, m), j] +
s[i, ifelse(j+1 ≤ n, j+1, 1)] + s[i, ifelse(j-1 ≥ 1, j-1, n)]
end
function ising2d_sweep!(s, β, niters)
m, n = size(s)
prob = [exp(-2*β*k) for k in -4:4]
for iter in 1:div(niters, m*n)
for j in 1:n, i in 1:m
s₁ = s[i,j]
k = s₁ * ising2d_sum_of_adjacent_spins(s, m, n, i, j)
s[i,j] = ifelse(rand() < prob[k+5], -s₁, s₁)
end
end
end
s = rand_ising2d(100)
@time ising2d_sweep!(s, β_crit, 10^9)
s
8.391483 seconds (218.71 k allocations: 11.477 MiB, 0.06% gc time)
100×100 Array{Int64,2}: -1 -1 1 -1 -1 -1 -1 1 1 1 … 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 -1 -1 1 -1 -1 -1 -1 -1 1 -1 -1 1 1 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 1 1 -1 -1 -1 1 1 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 1 1 -1 -1 -1 1 -1 1 1 -1 1 1 1 1 1 -1 -1 -1 -1 1 1 -1 1 -1 1 -1 -1 1 -1 … 1 1 1 1 1 -1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 1 -1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 … 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 -1 -1 -1 -1 -1 -1 -1 ⋮ ⋮ ⋱ ⋮ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 … 1 1 1 -1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 -1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 -1 -1 -1 -1 1 1 -1 -1 1 1 1 1 1 1 1 … 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 1 1 -1 -1 -1 -1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 1 1 -1 -1 -1 -1 -1 -1
# 2D Ising Metropolis (Gen Kuroki, 2017-01-09)
using PyPlot
const β_crit = log(1+sqrt(2))/2
rand_ising2d(n) = rand([-1, 1], n, n)
function ising2d_sum_of_adjacent_spins(s, m, n, i, j)
s[ifelse(i+1 ≤ m, i+1, 1), j] + s[ifelse(i-1 ≥ 1, i-1, m), j] +
s[i, ifelse(j+1 ≤ n, j+1, 1)] + s[i, ifelse(j-1 ≥ 1, j-1, n)]
end
function ising2d_sweep!(s, β, niters)
m, n = size(s)
prob = [exp(-2*β*k) for k in -4:4]
for iter in 1:div(niters, m*n)
for j in 1:n, i in 1:m
s₁ = s[i,j]
k = s₁ * ising2d_sum_of_adjacent_spins(s, m, n, i, j)
s[i,j] = ifelse(rand() < prob[k+5], -s₁, s₁)
end
end
end
function plot_ising2d(fig, s)
m, n = size(s)
pcolormesh(0:m, 0:n, s, vmin=-2.0, vmax=1.0, cmap="gist_earth")
fig[:set_aspect]("equal")
end
figure(figsize=(4,2))
s = rand_ising2d(100)
fig = subplot(121); plot_ising2d(fig, s); title("t = 0")
@time ising2d_sweep!(s, β_crit, 10^9); sleep(0.1)
fig = subplot(122); plot_ising2d(fig, s); title("t = 10\$^9\$")
tight_layout()
7.560066 seconds (112.12 k allocations: 5.614 MiB)
┌ Warning: `getindex(o::PyObject, s::Symbol)` is deprecated in favor of dot overloading (`getproperty`) so elements should now be accessed as e.g. `o.s` instead of `o[:s]`. │ caller = plot_ising2d(::PyCall.PyObject, ::Array{Int64,2}) at In[3]:29 └ @ Main .\In[3]:29
# 2D Ising Metropolis (Gen Kuroki, 2017-01-09)
# Pythonのmatplotlib.pyplotでプロットする
using PyPlot
# 2D Isingの臨界逆温度 (2D Isingはexactに解けている)
const β_crit = log(1+sqrt(2))/2
# ランダムに±1が配置されたn×nの2次元配列を返す函数
# この函数は 2D Ising model のシミュレーションの初期値として使われる
rand_ising2d(n) = rand([-1, 1], n, n)
# m×nの2次元配列 s を周期境界条件で考える.
# s[i,j] の値をスピンと呼ぶ.
# 次は (i,j) 成分の上下左右の4つのスピンの和を返す函数.
function ising2d_sum_of_adjacent_spins(s, m, n, i, j)
s[ifelse(i+1 ≤ m, i+1, 1), j] + s[ifelse(i-1 ≥ 1, i-1, m), j] +
s[i, ifelse(j+1 ≤ n, j+1, 1)] + s[i, ifelse(j-1 ≥ 1, j-1, n)]
end
# ここで ifelse(a, b, c) は a が真のとき b を, そうでないとき c を返す函数である.
# Julia言語では配列のインデックスが1から始まることに注意!
# 以下はメトロポリス法で m×n の配列 s の全体を niters÷(mn) 回更新する函数.
# βは逆温度と呼ばれる正の実数.
function ising2d_sweep!(s, β, niters)
m, n = size(s) # s は m×n の配列であるとする (s の成分は±1のみであることを仮定)
prob = [exp(-2*β*k) for k in -4:4]
# ↑ 繰り返し使われる確率計算の値 exp(-2βk), k ∈ {-4,-3,...,3,4} を配列に格納
# Julia言語の配列のインデックスは1から始まるので prob[k+5] = exp(-2βk) になる.
for iter in 1:div(niters, m*n)
for j in 1:n, i in 1:m
s₁ = s[i,j]
k = s₁ * ising2d_sum_of_adjacent_spins(s, m, n, i, j)
s[i,j] = ifelse(rand() < prob[k+5], -s₁, s₁)
# ↑ k を s[i,j] と周囲4つのスピンの和の積とするとき
# s[i,j] を -1 倍する確率は prob[k+4] = exp(-2βk) になる.
# 正確には確率は min(1.0, prob[k+5]) であるのだが、
# 更新の判定条件を rand() < prob[k+4] としているので同じことである.
end
end
end
# PyPlot を使って2次元配列 s をプロットするための函数
function plot_ising2d(fig, s)
m, n = size(s)
pcolormesh(0:m, 0:n, s, vmin=-2.0, vmax=1.0, cmap="gist_earth")
fig[:set_aspect]("equal")
end
figure(figsize=(4,2)) # プロット開始、画像サイズは横×縦=4inch×2inch
s = rand_ising2d(100) # 初期状態をランダムに設定
fig = subplot(121); plot_ising2d(fig, s); title("t = 0") # 初期状態をプロット
@time ising2d_sweep!(s, β_crit, 10^9); sleep(0.1) # 10億点を更新
fig = subplot(122); plot_ising2d(fig, s); title("t = 10\$^9\$") # 更新結果をプロット
tight_layout() # 二枚のプロットの端が重ならいように調節する
7.723458 seconds (112.11 k allocations: 5.614 MiB, 0.79% gc time)
┌ Warning: `getindex(o::PyObject, s::Symbol)` is deprecated in favor of dot overloading (`getproperty`) so elements should now be accessed as e.g. `o.s` instead of `o[:s]`. │ caller = plot_ising2d(::PyCall.PyObject, ::Array{Int64,2}) at In[4]:47 └ @ Main .\In[4]:47