Skip to content

Instantly share code, notes, and snippets.

@enakai00
Last active January 4, 2023 08:35
Show Gist options
  • Save enakai00/9aace6551be63848e9ce to your computer and use it in GitHub Desktop.
Save enakai00/9aace6551be63848e9ce to your computer and use it in GitHub Desktop.
Gibbs Sampling Emulation of 2D Ising Model
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pandas import Series, DataFrame
from numpy.random import randint, randn, rand
Size = 50
J = 1
H = 0.0
Temp = 0
def spin_direction(field, x, y):
energy = H
for dx, dy in [(-1,0), (1,0), (0,-1), (0,1)]:
# Cyclic boundary condition
if x + dx < 0: dx += Size
if y + dy < 0: dy += Size
if x + dx >= Size: dx -= Size
if y + dy >= Size: dy -= Size
energy += J * field.ix[x + dx, y + dy]
if Temp == 0:
p = (np.sign(energy) + 1) * 0.5
else:
p = 1/(1+np.exp(-2*(1/Temp)*energy))
if rand() <= p:
spin = 1
else:
spin = -1
return spin
def run_gibbs_sampling(field, iternum=5):
for _ in range(iternum):
lattice = DataFrame([(y,x) for x in xrange(Size) for y in xrange(Size)])
lattice.reindex(np.random.permutation(lattice.index))
for x, y in lattice.values:
field.ix[x, y] = spin_direction(field, x, y)
if __name__ == '__main__':
fig = plt.figure()
field = DataFrame(randint(2,size=(Size,Size))*2-1)
temps = [0.,.5,1.,1.5,2.,2.5,5.0,10.0][::-1]
for i in range(1,9):
Temp = temps[i-1]
run_gibbs_sampling(field)
ax = fig.add_subplot(2,4,i)
ax.set_title("T = %2.1f\nH = %1.1f" % (Temp, H))
axim = ax.imshow(field.values, vmin=-1, vmax=1,
cmap=plt.cm.gray_r, interpolation='nearest')
fig.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment