Juliaで学ぶ古典モンテカルロシミュレーション

おまけ:書き方による速度の違い

せっかくなので、一番コンパクトな書き方をしてみよう。そして、少しづつ書き方を変えてみる。速度差は出るだろうか。
まず、そのまま書き下してみる。 計算環境はJulia Boxを使用
https://www.juliabox.com
し、計算機環境は

In [1]:
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)

となっている。
二次元イジング模型のメトロポリス法によるモンテカルロシミュレーションのコードは、

In [2]:
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
Out[2]:
モンテカルロ法! (generic function with 1 method)

となり、このくらい短く書くことができる。

In [3]:
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)
Out[3]:
0.9920440608

次に、乱数を少し先に確保してみる。

In [4]:
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
Out[4]:
スイープ! (generic function with 1 method)
In [5]:
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)
Out[5]:
0.988445197

少し遅くなった? 念のためもう一回実行する。

In [6]:
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)
Out[6]:
0.988445197

やはり遅くなっている。配列を確保したのが原因だろうか?

乱数に関しては元に戻して、次に、ifelse関数を三項演算子「a ? b : c 」に変えてみる

In [7]:
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
Out[7]:
スイープ! (generic function with 1 method)
In [8]:
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)
Out[8]:
0.9920440608

遅い気がする。念のためもう一度。

In [9]:
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)
Out[9]:
0.9920440608

うーん。これは遅い。つまり、現時点では

In [10]:
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
Out[10]:
スイープ! (generic function with 1 method)
In [11]:
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)
Out[11]:
0.9920440608

が一番速い。
次に、ifelse関数を余りを計算するmod1(ix-1,Lx)で置き換えてみよう。この関数は

In [12]:
ix = 1
println(mod1(ix-1,Lx))
ix = Lx
println(mod1(ix+1,Lx))
100
1

と周期境界条件を表現できる。これを使ってみる。

In [13]:
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
Out[13]:
スイープ! (generic function with 1 method)
In [14]:
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)
Out[14]:
0.9920440608

とても遅くなった。ifelse関数が一番速い。

次に、コードが少し汚くなるが、ifelse関数も除去してみよう。

In [15]:
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
Out[15]:
スイープ! (generic function with 1 method)
In [16]:
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)
Out[16]:
0.9990793904

すごい汚いコードになったが速くはなった。どうせならさらにいってみよう。mzの部分はもう少し高速化できる。

In [17]:
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
Out[17]:
スイープ! (generic function with 1 method)
In [18]:
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)
Out[18]:
0.9990793904

これは余り速くならなかった。しかし、あえてifelse関数を使う必要はないので、この変更は採用しよう。
結局、見やすさと速さを両立しているのは、

In [2]:
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
Out[2]:
モンテカルロ法! (generic function with 1 method)
In [20]:
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)
Out[20]:
0.9920440608

この形が良さそうだ。

In [1]:
using Plots
pyplot()
Out[1]:
Plots.PyPlotBackend()
In [ ]:
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
In [ ]:
plot(time,Mztdep)