improving python code in Monte Carlo simulation

Monte Carlo simulation

At first I get rid of the lists, afterwards I used a just in time compiler (numba). Without compilation a get 196s (your version), with compilation I get 0.44s. So there is a improvement of a factor 435 by using a jit-compiler and a few easy optimizations here.

A main advantage is also, that the GIL (global interpreter lock) is released here also. If the code is processor limited and not limited by the memory bandwith, the random numbers can be calculated in another thread while running the simulation in the other (more than one core can be used). This could also improve performance a bit further and would work as follows:

1. Calculate the first chunk of random numbers (Small enough that the whole problem fits easily in the processor cache (at least L3 cache).
2. Start an iteration with that random numbers. While running the iteration another chunk of random numbers is calculated.
3. Repeat (2) until you are done.

Code

import numba as nb
import numpy as np

def calc (m,j,e,r,dt,b,rgrid,mgrid):
    M=m*m
    N = M * r
    i=0
    while i < j:
        # 1. platz aussuchen
        x = np.random.randint(M)  # wähle zufälligen platz aus

        if rgrid[x] != 0:
            i += dt / N

            ########
            # 2. teilchen aussuchen
            #li1 = np.arange(-abs(rgrid[x]), abs(rgrid[x]) + 1, 2)
            #li2 = np.arange(0, abs(rgrid[x]) + 1)

            #li3 = zip(li1, li2)  # list1 und list2 als tupel in list3
            #results = [item[1] for item in li3 if item[0] == mgrid[x]]  # gebe 2. element von tupel aus für passende magnetisierung
            #num = int(''.join(map(str, results)))  # wandle listeneintrag in int um
            #######

            # This should be equivalent
            jj=0 #li2 starts with 0 and has a increment of 1

            for ii in range(-abs(rgrid[x]),abs(rgrid[x])+1, 2):
                if (ii==mgrid[x]):
                    num=jj
                    break

                jj+=1

            spin = 1.0 if np.random.random() < num / rgrid[x] else -1.0

            # 3. ereignis aussuchen
            p = np.random.random()
            p1 = np.exp(- spin * b * mgrid[x] / rgrid[x]) * dt  # flip
            p2 = dt  # hoch
            p3 = dt  # runter
            p4 = (1 - spin * e) * dt  # links
            p5 = (1 + spin * e) * dt  # rechts
            p6 = 1 - (4 + p1) * dt  # nichts


            if p < p6:
                continue
            elif p < p6 + p1:
                #flip()
                mgrid[x] -= 2 * spin 
                continue
            elif p < p6 + p1 + p2:
                #up()
                y = x - m
                if y < 0:  # periodische randbedingung hoch
                    y += m * m
                rgrid[x] -= 1  # [zeile, spalte] masse -1
                rgrid[y] += 1  # [zeile, spalte] masse +1
                mgrid[x] -= spin  # [zeile, spalte] spinänderung alter platz
                mgrid[y] += spin  # [zeile, spalte] spinänderung neuer platz
                continue
            elif p < p6 + p1 + p2:
                #down()
                y = x + m
                if y >= m * m:  # periodische randbedingung unten
                    y -= m * m
                rgrid[x] -= 1
                rgrid[y] += 1
                mgrid[x] -= spin
                mgrid[y] += spin
                continue
            elif p < p6 + p2 + p3:
                #left()
                y = x - 1
                if (y // m) < (x // m):  # periodische randbedingung links
                    y += m
                rgrid[x] -= 1
                rgrid[y] += 1
                mgrid[x] -= spin
                mgrid[y] += spin
                continue
            else:
                #right()
                y = x + 1
                if (y // m) > (x // m):  # periodische randbedingung rechts
                    y -= m
                rgrid[x] -= 1
                rgrid[y] += 1
                mgrid[x] -= spin
                mgrid[y] += spin
                continue
    return (mgrid,rgrid)

And now the main function for testing:

def main():
    m = 10  # zeilen, spalten
    j = 1000 # finale zeit
    r = 3  # platzdichte
    b = 1.6  # beta
    e = 0.9  # epsilon

    M = m * m  # platzanzahl
    N = M * r  # teilchenanzahl
    dt = 1 / (4 * np.exp(b))  # delta-t
    i = 0

    rgrid = r * np.ones((m, m),dtype=np.int) #don't convert the array build it up with the right datatype  # dichte-matrix, rho = n(+) + n(-)
    magrange = np.arange(-r, r + 1, 2)  # mögliche magnetisierungen [a, b), schrittweite
    mgrid = np.random.choice(magrange, (m, m))  # magnetisierungs-matrix m = n(+) - (n-)

    #Compile the function
    nb_calc = nb.njit(nb.types.Tuple((nb.int32[:], nb.int32[:]))(nb.int32, nb.int32,nb.float32,nb.int32,nb.float32,nb.float32,nb.int32[:], nb.int32[:]),nogil=True)(calc)

    Results=nb_calc(m,j,e,r,dt,b,rgrid.flatten(),mgrid.flatten())

    #Get the results
    mgrid_new=Results[0].reshape(mgrid.shape)
    rgrid_new=Results[1].reshape(rgrid.shape)

Edit
I have rewritten the code “2.Teilchen aussuchen” and reworked the code so that it works with scalar indices. This gives an additional speed up by a factor of 4.