Snippet:

import numpy
import math
import random
def init_isig(Nx, Ny, random=True):
if random:
return 2*numpy.random.randint(2, size=(Nx, Ny)) -1
else:
return numpy.ones(shape=(Nx,Ny))
def delta_energy(U, x, y, b, h):
(Nx, Ny) = U.shape
x_l = (x + Nx -1) % Nx
x_h = (x + Nx +1) % Nx
y_l = (y + Ny -1) % Ny
y_h = (y + Ny +1) % Ny
E = 2 * U[x,y] * (h + b * (U[x_l,y]+U[x_h,y]+U[x,y_l]+U[x,y_h]))
return E
def iseg_metropolis(U, b, h):
(Nx, Ny) = U.shape
x = random.randint(0,Nx-1)
y = random.randint(0,Ny-1)
dE = delta_energy(U, x, y, b, h)
if random.random() < math.exp(-dE):
U[x,y] = -U[x,y]
return U
if __name__ == "__main__":
U = init_isig(10, 10)
b = -1.0
h = 0
for _ in range(100*100):
U = iseg_metropolis(U, b, h)
print(U)