1

For some reason there is an increase of execution time around (10:1) (4min vs 40min) in the following code: Code 1

def E_site(lattice, i, j):
    N = len(lattice)
    Nb = lattice[(i+1)%N,j] + lattice[(i-1)%N,j] + lattice[i,(j+1)%N] + lattice[i, (j-1)%N]
    return -lattice[i, j]*Nb

def metropolis(lattice, T, Neq):
    for n in range(Neq):
        N = len(lattice)
        i = np.random.randint(0, N)
        j = np.random.randint(0, N)
        Ei = E_site(lattice, i, j)
        lattice[i, j] *= -1
        Ef = E_site(lattice, i, j)
        dE = Ef - Ei
        if dE < 0 or  np.random.rand() < np.exp(-dE/T):
            pass
        else:
            lattice[i,j] *= -1
    return lattice
def equilibrate(lattice, T, N_equilibration, show_stats):
    if show_stats == False:
        lattice = metropolis(lattice, T, N_equilibration)
        return lattice

def simulate_temperatures(Ti, Tmax, Tstep, N_eq, N_post):
    avg_m_list, avg_E_list = [], []
    T = Ti
    while(T <= Tmax):
        s = time.clock()
        if T <= 2:
            N_eq = int(1e3)
            lattice = init_lattice(20,0.5)
            eq_lattice = equilibrate(lattice, T, N_eq, False)
        else:
            lattice = init_lattice(20,0.5)
            eq_lattice = equilibrate(lattice, T, N_eq, False)
        E, m = [], []
        for i in range(N_post):
            lattice = metropolis(eq_lattice,T,N_eq)
            E.append(total_avg_E(lattice))
            m.append(calcMag(lattice))
        T += Tstep

Code 2

    def E_site(self, i, j):
        N = len(self.lattice)
        Nb = self.lattice[(i+1)%N,j] + self.lattice[(i-1)%N,j] + self.lattice[i,(j+1)%N] +     self.lattice[i, (j-1)%N]
        return -self.lattice[i, j]*Nb

    def alternative_metropolis(self, Neq, T):
        for n in range(Neq):
            N = len(self.lattice)
            i = np.random.randint(0, N)
            j = np.random.randint(0, N)
            Ei = self.E_site(i, j)
            self.lattice[i, j] *= -1
            Ef = self.E_site(i, j)
            dE = Ef - Ei
            if dE < 0 or  np.random.rand() < np.exp(-dE/T):
                pass
            else:
                self.lattice[i,j] *= -1


    def alternative_several_temperatures(self, Ti, Tf):
        Tlist = np.arange(Ti, Tf, 0.1)
        for T in Tlist:
            s = time.clock()
            for i in range(len(self.lattice)):
                for j in range(len(self.lattice)):
                    if random() < 0.49: 
                        self.lattice[i,j] = 1
                    else:
                        self.lattice[i,j] = -1
            self.alternative_metropolis(int(1e7), T)

I put all the functions that are called at some point to let you check they look the same, except Code 2 is inside a class. Then, in Code1 > simulate_temperatures > for loop with Npost there is a call to metropolis. The inside metropolis loop, will be done then Npost*Neq times. So if Npost = 1e4, Neq = 1e3 then the inside metropolis loop is done 1e7 times in total.

If you go to the method alternative_several_temperatures in Code 2 you see that the alternative_metropolis method has an argument 1e7 , that is the number of times the loop inside alternative_metropolis will be executed. Therefore, as all the functions look the same to me, how is it Code 2 runs at 40min whereas Code 1 runs at 4min with Npost = 1e4, Neq = 1e3 ? Did i do something wrong with the math?

*Actually Code 2 does more iterations due to equilibrate function calls.

bp41
  • 137
  • 13
Jay
  • 95
  • 8
  • 1
    Referring to the latter part, could you rephrase it and make it more readable ? – bp41 Jun 23 '20 at 19:19
  • People might offer suggestions, but really the best answer and the first line of attack for questions like this is always going to be to profile your code: https://stackoverflow.com/questions/582336/how-can-you-profile-a-python-script. – sgillen Jun 23 '20 at 19:32

1 Answers1

0

There's a lot going on here, so it's difficult to tell with certainty, but I assume you're calling this code thousands of times per second. If that's the case, then this is likely to be the culprit:

    i = np.random.randint(0, N)
    j = np.random.randint(0, N)

The phenomenon is known as "entropy pool exhaustion." Essentially, your computer generates entropy by accumulating entropy from a variety of sources (keyboard strokes, mouse movement, ambient temperature, et al.). The rate at which this entropy accumulates is rather slow. When you request additional random integers and the entropy pool has already been exhausted, the machine blocks, waiting for more entropy to become available.

To solve this you need to use a smaller subset of "true" random numbers and a pseudorandom number generator to fill in the gaps. You could, for example, create an "entropy provider" class which re-seeds its PRNG with a new random integer only every 10,000 calls.

Max
  • 913
  • 1
  • 7
  • 18
  • Usually if a random number generator is using an entropy pool it will be mentioned, otherwise a pseudorandom generator is assumed. I don't see any mention of it in the [`randint` documentation](https://numpy.org/devdocs/reference/random/generated/numpy.random.randint.html). That said, a good pseudorandom generator can be slower than you expect. – Mark Ransom Jun 23 '20 at 20:11