Ising model cython

cyising.pyx

# Reference: http://nbviewer.jupyter.org/gist/genkuroki/4fa46c68c56ee0f3b1a6fc8ec628b9d7
# Author MathSorcerer

from math import exp
import array
from cpython cimport array as carray
from random import choice, random
from libc.stdlib cimport rand, RAND_MAX
import time
cimport cython
import numpy as np

@cython.boundscheck(False)
cdef int ising2d_sum_of_adjacent_spins(int[:,:] s, int m, int n,int  i, int j):
    cdef:
        int i_bottom = i+1 if i+1 < m else 0
        int i_top = i-1 if i-1 >= 0 else m-1
        int j_right = j+1 if j+1 < n else 0
        int j_left = j-1 if j-1 >= 0 else n-1
    return s[i_bottom,j]+s[i_top,j]+s[i,j_right]+s[i,j_left]

@cython.boundscheck(False)
@cython.wraparound(False)
cpdef int[:,:] ising2d_sweep(int[:,:] s, double beta, int niters):
    cdef int m, n, s1, loop,num_iteration,k
    m,n= s.shape[0],s.shape[1]
    num_iteration=int(niters/(m*n))
    cdef carray.array prob = array.array('d',[exp(-2*beta*k) for k in [-4, -3, -2, -1, 0, 1, 2, 3, 4]])
    for loop in range(num_iteration):
        for i in range(m):
            for j in range(n):
                s1 = s[i,j]
                k = s1*ising2d_sum_of_adjacent_spins(s, m, n, i, j)
                s[i,j] = -s1 if (rand()+1.0)/(RAND_MAX+2.0) <  prob[k+4] else s1
    return s

cymain.py

import pyximport 
pyximport.install()
from math import log,sqrt
import numpy as np 
from cyising import ising2d_sweep
from random import choice
from copy import deepcopy
import time

from matplotlib import pyplot as plt 

def plot_result(s_begin,s_end):
    fig = plt.figure()
    ax1 = fig.add_subplot(121)
    ax2 = fig.add_subplot(122)
    ax1.imshow(s_begin, cmap='gray')
    ax2.imshow(s_end, cmap='gray')
    plt.show()

from copy import deepcopy
from math import log,sqrt 
import numpy as np
n=100
beta_critical = log(1+sqrt(2))/2

def main():
    rand_ising2d = np.array([[choice([-1, 1]) for j in range(n)] for i in range(n)]).astype(np.int32)
    s_begin = deepcopy(rand_ising2d)
    begin = time.time()
    s_end = ising2d_sweep(rand_ising2d, beta_critical, 1e9)
    end = time.time()
    print("Elapsed=", end-begin)
    plot_result(s_begin,s_end)

if __name__ == '__main__':
    main()

gistfile1.txt

import pyximport 
pyximport.install()
from math import log,sqrt
import numpy as np 
from cyising import ising2d_sweep
from random import choice
from copy import deepcopy
import time

from matplotlib import pyplot as plt 

def plot_result(s_begin,s_end):
    fig = plt.figure()
    ax1 = fig.add_subplot(121)
    ax2 = fig.add_subplot(122)
    ax1.imshow(s_begin, cmap='gray')
    ax2.imshow(s_end, cmap='gray')
    plt.show()

from copy import deepcopy
from math import log,sqrt 
import numpy as np
n=100
beta_critical = log(1+sqrt(2))/2

def main():
    rand_ising2d = np.array([[choice([-1, 1]) for j in range(n)] for i in range(n)]).astype(np.int32)
    s_begin = deepcopy(rand_ising2d)
    begin = time.time()
    s_end = ising2d_sweep(rand_ising2d, beta_critical, 1e9)
    end = time.time()
    print("Elapsed=", end-begin)
    plot_result(s_begin,s_end)

if __name__ == '__main__':
    main()

License is MIT