03_IsingMCMC.diff

03_IsingMCMC.diff

--- ../03_IsingMCMC.jl-master/IsingMCMC02a.jl	2020-07-05 09:45:56.000000000 +0900
+++ IsingMCMC02a.jl	2020-07-05 19:34:28.617972200 +0900
@@ -1,6 +1,8 @@
 module IsingMCMC
 ##export E, MCMC
 
+using Random: default_rng
+
 function genFunctor( getJ, getB)
     getJ = getJ
     getB = getB
@@ -10,7 +12,8 @@
         Et = 0
         for i in 1:N
             J = Main.getNeighbors(i)
-            Et += -0.5*S[i]*sum( [getJ(i,j) * S[j] for j in J] )
+            #Et += -0.5*S[i]*sum(getJ(i,j) * S[j] for j in J)
+            Et += -0.5*S[i]*sum(S[j] for j in J)
             Et += -    S[i]*getB(i)
         end
         Et
@@ -20,7 +23,8 @@
     function dE(S, k, N)
         dEt = 0
         J = Main.getNeighbors(k)
-        dEt +=   S[k]*sum( [(getJ(k,j) + getJ(j,k)) * S[j] for j in J] )
+        #dEt +=   S[k]*sum((getJ(k,j) + getJ(j,k)) * S[j] for j in J)
+        dEt += 2*S[k]*sum(S[j] for j in J)
         dEt += 2*S[k]*getB(k)
         dEt
     end
@@ -36,26 +40,23 @@
     # MCMC
     function MCMC( T, N, trial)
         # initialize
-        sim = []
-        S = ones(N)
+        S = ones(Int8, N)
+        sim = typeof(S)[]
         for i in 1:N
             if rand() < 0.5
                 flipx!(S,i) # random flip at first
             end
         end
-        # Gibbs sampling position
-        gpos = [rand(1:N) for i in 1:trial]
-        # random value in [0,1]
-        rn = [rand() for i in 1:trial]
-
+        
         # MCMC trial
+        rng = default_rng()
         for t in 1:trial
-            k = gpos[t]
+            k = rand(rng, 1:N) 
             de = dE(S, k, N)
-            if rn[t] < PrE(de, T) # MH criteria 
+            if rand(rng) < PrE(de, T) # MH criteria 
                 flipx!(S, k) # change k
             end
-            push!(sim, copy(S))
+            iszero(t % 1000) && push!(sim, copy(S))
         end
 
         sim
--- ../03_IsingMCMC.jl-master/IsingSimMain02a.jl	2020-07-05 09:45:56.000000000 +0900
+++ IsingSimMain02a.jl	2020-07-05 18:40:52.711198600 +0900
@@ -13,13 +13,14 @@
 ###################################################################################################
 # view functions
 function drawEnergy(sim, trial, T, d=d, N=N)
-    xs = [t for t in 1:1000:trial]
+    #xs = [t for t in 1:1000:trial]
+    xs = 1:(trial ÷ 1000)
     ys = [E(sim[t],N) for t in xs]
     p = Plots.plot(xs, ys, label="Energy history with T=$(T)")
 end
 
 function drawSim(sim, f, d=d, N=N)
-    snap = sim[f]
+    snap = sim[f ÷ 1000]
     
     p=Plots.plot([],[],legend=false, size=(500,500))
     

https://qiita.com/oki_mebarun/items/25cf371601e8075458f9

で紹介されているコードにおける無駄なメモリ消費を減らすパッチ.

1. 函数 E と dE の中で getJ を使っている点が無駄だと考えられる. iとjが繋がっていれば get(i, j) は 1 になる.

2. 途中の10万個の状態をすべて記録に残している点も無駄だと考えらえる.

3. 状態Sの成分をInt8型にした.

という感じで実行すると、オリジナルでは大量にメモリを消費してガベージコレクションがひどいことになっている：

変更後はこうなる：