4
\$\begingroup\$

I am new to this community; I have tried my best to respect the policy of the community. I have written the Monte Carlo metropolis algorithm for the ising model. I want to optimize the code. I have tried my best. I want to optimize it further. The following is the code:

(I have used tricks like finding exponential only once, careful generation of random number, etc.)

import numpy as np
import time
import random

def monteCarlo(N,state,Energy,Mag,Beta,sweeps):
    if sweeps > 10000:
        print("Warning:Number of sweeps exceeded 10000\n\tsetting number of sweeps to 10000")
        sweeps = 10000
    start_time = time.time()
    expBeta = np.exp(-Beta*np.arange(0,9))
    E = np.zeros(sweeps)
    M = np.zeros(sweeps)
    for t in range(sweeps):
        for tt in range(N*N):
            a = random.randint(0, N-1)
            b = random.randint(0, N-1)
            s = state[a,b]
            delta_E = 2*s*(state[(a+1)%N,b] + state[a,(b+1)%N] + state[(a-1)%N,b] + state[a,(b-1)%N])
            if delta_E < 0:
                s *= -1
                Energy += delta_E
                Mag += 2*s
            elif random.random() < expBeta[delta_E]:
                s *= -1
                Energy += delta_E
                Mag += 2*s
            state[a, b] = s
        E[t] = Energy
        M[t] = Mag
    print("%d monte carlo sweeps completed in %d seconds" %(sweeps,time.time()-start_time))
    return E,M #returning list of Energy and Magnetization set

#####lattice config#####
"""N is lattice size
nt is number of Temperature points
sweeps are number of mc steps per spin"""
print("Starting Ising Model Simulation")
N = int(input("Enter lattice size : "))
startTime = time.time()
nt = 10
N2 = N*N
sweeps = 10000 #mc steps per spin
"""we will plot the following wrt temperature, T"""
T = np.linspace(2, 3, nt) #you can set temperature range
"""preparing lattice with all spins up"""
state = np.ones((N,N),dtype="int")
Energy = -N2
Mag = N2
#temperature loop
#for k in tqdm_gui(range(nt)):
for k in range(nt):
    temp = T[k]
    Beta=1/temp
    print("____________________________________\nTemperature is %0.2f, time is %d" %(temp,time.time()-startTime))
    E, M = monteCarlo(N,state,Energy,Mag,Beta,sweeps) #list of Energy and Magnetization
    Energy = E[-1]
    Mag = M[-1]

#further code is finding magnetization, autocorrelation, specific heat, autocorrelation, etc.
```
\$\endgroup\$
2
\$\begingroup\$

Naming variables

The usual rule of variable names in snake_case applies, i.e. energyFunctional would become energy_functional. Class names on the other hand should be written in CamelCase. I don't mind using single capital letters for matrices.

Why use a,b for discrete indices? I would use one of i,j,k,l,n,m,p,q,r.

Use descriptive names: e.g. energies instead of E.

Merging conditions

Instead of

if delta_E < 0:
    s *= -1
    Energy += delta_E
    Mag += 2*s
elif random.random() < expBeta[delta_E]:
    s *= -1
    Energy += delta_E
    Mag += 2*s

you could simply write

if delta_E < 0 or random.random() < expBeta[delta_E]:
    s *= -1
    Energy += delta_E
    Mag += 2*s

which is easier to read.

String formatting

Use the sweet f-strings.

sweep_time = int(time.time() - start_time)
print(f"{sweeps} monte carlo sweeps completed in {sweep_time} seconds")

Logging warnings

Consider using the logging library. I'd write warnings to stderr, but it's up to you, see.

import sys

print("Warning: Number of sweeps exceeded 10000", file=sys.stderr)
print("         setting number of sweeps to 10000", file=sys.stderr)

Truncating it to a single line allows for easier parsing later.

print("Warning: Number of sweeps truncated to 10000.", file=sys.stderr)

Refactorisation

If performance wasn't the primary goal, I'd introduce a few new functions. I'd separate the timing part anyway.

start_time = time.time()
energies, mags = monte_carlo(n, state, energy, mag, beta, sweeps)
elapsed_seconds = int(time.time() - start_time)
print(f"{sweeps} monte carlo sweeps completed in {elapsed_seconds} seconds")

monte_carlo

Applying these ideas, the monteCarlo function becomes the following.

def monte_carlo(n, state, energy, mag, beta, sweeps):
    if sweeps > 10000:
        print("Warning: Number of sweeps truncated to 10000.", file=sys.stderr)
        sweeps = 10000

    exp_betas = np.exp(-beta*np.arange(0,9))

    energies = np.zeros(sweeps)
    mags = np.zeros(sweeps)
    for t in range(sweeps):
        for tt in range(n*n):
            j = random.randint(0, n-1)
            k = random.randint(0, n-1)
            s = state[j,k]

            neighbour_sum = (state[(j-1)%n, k] +
                             state[j, (k-1)%n] + state[j, (k+1)%n] +
                             state[(j+1)%n, k])
            energy_diff = 2*s*neighbour_sum

            if energy_diff < 0 or random.random() < exp_betas[energy_diff]:
                s *= -1
                energy += energy_diff
                mag += 2*s

            state[j, k] = s

        energies[t], mags[t] = energy, mag

    return energies, mags
| improve this answer | |
\$\endgroup\$
  • \$\begingroup\$ Albeit, numpy random is better for scientific and statistical analysis purposes. (stackoverflow.com/a/7030595/8081272) @Gucio I have gone through your answer and you're correct. The only problem I am facing right now is, if I run my earlier code (the one in question), the average magnetisation is zero above temp >2.269. This ain't case when I use numpy's random. What I will do is, I will try @Andrew 's code with MT as generator. \$\endgroup\$ – Kartik Chhajed May 15 at 9:04
  • \$\begingroup\$ I have updated my question, pointing the problem faced with your suggestions wrt my original code. Though, your suggestion is very helpful for beginners. For example, as a newb coder, I earlier did not know when to use snake_case and when to use camelCase. \$\endgroup\$ – Kartik Chhajed May 15 at 11:31
  • \$\begingroup\$ @KartikChhajed Can you post the exact code used in generating your plots? Pastebin or something would work. \$\endgroup\$ – Andrew May 15 at 11:48
  • \$\begingroup\$ I have updated my question again with code. Sorry, I have not incorporated your suggestions in my code: like using snake_case. I was avoiding posting this complicated stuff earlier. Thank You for your advise given more first. \$\endgroup\$ – Kartik Chhajed May 15 at 12:07
  • \$\begingroup\$ @KartikChhajed No problem. Generally it is best not to edit your original question as the answers arrive. It is also good practice to not accept answers too quickly. It seems that I have overlooked something regarding the batch-generation of these pseudo-random numbers. I've removed the relevant parts from my answer. I'll update it if I think of something useful. \$\endgroup\$ – Andrew May 15 at 13:30
4
\$\begingroup\$

I cannot yet comment due to my reputation, so I'll write this as an answer to your comment on Andrew's answer, and delete if someone comments this information or Andrew updates his answer.

Saying that numpy's random is not a good random number generator does not seem right to me. From numpy's reference :

By default, Generator uses bits provided by PCG64 which has better statistical properties than the legacy MT19937 used in RandomState.

So to me it seems that :

  1. numpy uses the PCG64 random generator which, according to numpy, has better statistical properties than legacy MT19937
  2. numpy used to use MT19937
  3. you can still chose to use the MT19937 random number generator
from numpy.random import Generator, MT19937
rg = Generator(MT19937(12345))
rg.random()

Maybe there is something I'm missing, and maybe it is in this part of your comment

It produces 53-bit precision floats and has a period of 2**19937-1

If so, I'd be interested to know how numpy's random would still be flawed for a Monte-Carlo analysis.

| improve this answer | |
\$\endgroup\$
  • \$\begingroup\$ I have tried importing from numpy.random import Generator. But, it shows import error: ImportError: cannot import name 'Generator' from 'numpy.random' (/usr/lib/python3/dist-packages/numpy/random/__init__.py). \$\endgroup\$ – Kartik Chhajed May 15 at 9:30
  • \$\begingroup\$ Generator seems to have appeared in numpy v1.17.0, so see if your numpy version is up to date. I can run the snippet on numpy v1.18.1. \$\endgroup\$ – Gucio May 15 at 9:37
  • \$\begingroup\$ An answer that elaborates on random number generation is always welcome. No need to delete yours. \$\endgroup\$ – Andrew May 15 at 10:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.