How can i speed up the MC simulation of Ising model?

Question

I am new here.
I programmed a Monte Carlo simultaion of 2D Ising model use Python, and find it is very slow. and i don't know how to speed it up.
I post the code below(only the test of trail_move function) and used line_profiler to count the running time of each line, i noticed that the random() take me almost 60% of the Time, I dont know if it can be solved.

import time, random
import numpy
import math


L = 10
lattice_num = L**2
MCSs = 10000

T=2.1

S_lattices = [[0 for col in range(L)] for row in range(L)]
for i in range(L):
    S_lattices[i][:]=random.choices([-1,1],weights=[1,1],k=L)
Exp_value = [1, math.exp(-4/T) , math.exp(-8/T)]

def trail_move():
    index_x = numpy.random.randint(0, L)
    index_y = numpy.random.randint(0, L)

    # delta = 'n' - 'o'
    delta_M = -2*S_lattices[index_x][index_y]
    S_adjacent = S_lattices[index_x][(index_y-1)%L] + S_lattices[index_x][(index_y+1)%L] + S_lattices[(index_x+1)%L][index_y] + S_lattices[(index_x-1)%L][index_y]
    delta_E = -delta_M * S_adjacent #delta_E = [-(-S) * S_sum-(-S*S_sum)]
    if delta_E < 0:
        S_lattices[index_x][index_y] = -S_lattices[index_x][index_y]
    else:
        random_num = numpy.random.random()
        Prob = Exp_value[int(delta_E//4)] # math.exp(-delta_E/T)
        if random_num <= Prob:
            S_lattices[index_x][index_y] = -S_lattices[index_x][index_y]

def init_move(sweep_time=MCSs):
    for i in range(sweep_time):
        for j in range(L**2):
            trail_move()

from line_profiler import LineProfiler
lprofiler = LineProfiler(init_move, trail_move)
lprofiler.run('init_move()')
lprofiler.print_stats()

And the Running time of each line is below

Timer unit: 1e-07 s

Total time: 9.56945 s
File: D:\my_research\Homework\MC_Ising\5optimization.py
Function: trail_move at line 27

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    27                                           def trail_move():
    28   1000000   28009338.0     28.0     29.3      index_x = numpy.random.randint(0, L)
    29   1000000   27092754.0     27.1     28.3      index_y = numpy.random.randint(0, L)
    30                                           
    31                                               # delta = 'n' - 'o'
    32   1000000    5287582.0      5.3      5.5      delta_M = -2*S_lattices[index_x][index_y]
    33   1000000   10942671.0     10.9     11.4      S_adjacent = S_lattices[index_x][(index_y-1)%L] + S_lattices[index_x][(index_y+1)%L] + S_lattices[(index_x+1)%L][index_y] + S_lattices[(index_x-1)%L][index_y]
    34   1000000    3764977.0      3.8      3.9      delta_E = -delta_M * S_adjacent #delta_E = [-(-S) * S_sum-(-S*S_sum)]
    35   1000000    3255163.0      3.3      3.4      if delta_E < 0:
    36     38443     205736.0      5.4      0.2          S_lattices[index_x][index_y] = -S_lattices[index_x][index_y]
    37                                               else:
    38    961557    8242272.0      8.6      8.6          random_num = numpy.random.random()
    39    961557    4844328.0      5.0      5.1          Prob = Exp_value[int(delta_E//4)] # math.exp(-delta_E/T)
    40    961557    3693935.0      3.8      3.9          if random_num <= Prob:
    41     67664     355778.0      5.3      0.4              S_lattices[index_x][index_y] = -S_lattices[index_x][index_y]

Total time: 13.0176 s
File: D:\my_research\Homework\MC_Ising\5optimization.py
Function: init_move at line 43

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    43                                           def init_move(sweep_time=MCSs):
    44     10001      28754.0      2.9      0.0      for i in range(sweep_time):
    45   1010000    2728568.0      2.7      2.1          for j in range(L**2):
    46   1000000  127419066.0    127.4     97.9              trail_move()

[Finished in 14.0s]

And I also find a course(here), post a problem to Speeding Up Your Code, give a Hint that I should work with many of the spins at once, but I still don't know how to deal with it.

a screenshot of the Hint in that website

Need help!!!
Q: How can i speed up my program? (Use the python program only)
q1: how can i work with many of spins at once?
q2: how can i deal with random function with so slow speed?

user_10429349 · Answer

I would use a numba decorator there. In my test it made the code faster by a factor ~30

import numpy as np
import math
import numba as nb

L = 10
lattice_num = L**2
MCSs = 10000

T=2.1

Exp_value = np.array([1, np.exp(-4/T) , np.exp(-8/T)])
S_lattices = np.zeros((L, L))
for i in range(L):
    S_lattices[i, :] = random.choices([-1,1], weights=[1,1], k=L)


@nb.njit
def init_move(S_lattices, sweep_time=MCSs):
    for i in range(sweep_time):
        for j in range(L**2):
            index_x = np.random.randint(0, L)
            index_y = np.random.randint(0, L)

            # delta = 'n' - 'o'
            delta_M = -2*S_lattices[index_x, index_y]
            S_adjacent = S_lattices[index_x, (index_y-1)%L] + S_lattices[index_x, (index_y+1)%L] + S_lattices[(index_x+1)%L, index_y] + S_lattices[(index_x-1)%L, index_y]
            delta_E = -delta_M * S_adjacent #delta_E = [-(-S) * S_sum-(-S*S_sum)]
            if delta_E < 0:
                S_lattices[index_x, index_y] = -S_lattices[index_x, index_y]
            else:
                random_num = np.random.random()
                Prob = Exp_value[int(delta_E//4)] # math.exp(-delta_E/T)
                if random_num <= Prob:
                    S_lattices[index_x, index_y] = -S_lattices[index_x, index_y]

init_move(S_lattices)

user_5140209 · Answer

I think working with many spins at once means either trying to update many spins at once, making sure that they are far apart enough that can be updated independently (like on a chessboard pattern) or using cluster moves such the Wolff algorithm. In both cases you need to re-think your code almost from scratch.

Having said that, you can still improve the performances of your code if you look into it a bit more deeply. It turns out that numpy is not very good at generating a small number of random integer numbers, and that in python in general it's more efficient to generate floating point random numbers and then casting them to integers. In addition, your code has a few redundant statements that can be joined together to further improve performances.

On my machine your code with the following trail_move function runs 2.5 times faster than the original one, and the time taken by trail_move itself decreases by a factor of roughly 3.5.

def trail_move():
    index_x = int(random.random() * L)
    index_y = int(random.random() * L)

    S_adjacent = S_lattices[index_x][(index_y-1)%L] + S_lattices[index_x][(index_y+1)%L] + S_lattices[(index_x+1)%L][index_y] + S_lattices[(index_x-1)%L][index_y]
    delta_E = 2*S_lattices[index_x][index_y] * S_adjacent
    if delta_E < 0 or random.random() <= Exp_value[int(delta_E//4)]:
        S_lattices[index_x][index_y] = -S_lattices[index_x][index_y]

How can i speed up the MC simulation of Ising model?

2 Answers