Hi, I have been writing a program that simulates the behaviour of a system of hard disks using the monte carlo method. I ran the program upto around 2 million steps and it's really slow(takes around 3 1/2 to 4 hours).

I first wrote it in vpython as I'm a complete beginner and I needed to see where I was going wrong!

However I'm trying to move away from the visual side now as that should help speed things up. I tested psyco over 5,000 steps and found that it was actually 4 seconds slower than without it!

I was just wondering if anyone could give me any pointers as to things that could obviously be optimised.

Here's the full code, Ignore the PairCorrelationFunction function as I didn't write that part.

``````from visual import *	#

from math import *	#IMPORT LIBRARIES

from random import *	#
from numpy import array

import time

start = time.time()

def PairCorrelationFunction_2D(x,y,S,rMax,dr):
"""Compute the two-dimensional pair correlation function, also known
as the radial distribution function, for a set of circular particles
contained in a square region of a plane.  This simple function finds
reference particles such that a circle of radius rMax drawn around the
particle will fit entirely within the square, eliminating the need to
compensate for edge effects.  If no such particles exist, an error is
returned. Try a smaller rMax...or write some code to handle edge effects! ;)

Arguments:
x               an array of x positions of centers of particles
y               an array of y positions of centers of particles
S               length of each side of the square region of the plane
rMax            outer diameter of largest annulus
dr              increment for increasing radius of annulus

Returns a tuple: (g, radii, interior_x, interior_y)
g(r)            a numpy array containing the correlation function g(r)
radii           a numpy array containing the radii of the
annuli used to compute g(r)
interior_x      x coordinates of reference particles
interior_y      y coordinates of reference particles
"""
from numpy import zeros, sqrt, where, pi, average, arange, histogram
# Number of particles in ring/area of ring/number of reference particles/number density
# area of ring = pi*(r_outer**2 - r_inner**2)

# Find particles which are close enough to the box center that a circle of radius
# rMax will not cross any edge of the box

bools1 = x>1.1*rMax
bools2 = x<(S-1.1*rMax)
bools3 = y>rMax*1.1
bools4 = y<(S-rMax*1.1)
interior_indices, = where(bools1*bools2*bools3*bools4)
num_interior_particles = len(interior_indices)

if num_interior_particles < 1:
raise  RuntimeError ("No particles found for which a circle of radius rMax\
will lie entirely within a square of side length S.  Decrease rMax\
or increase the size of the square.")

edges = arange(0., rMax+1.1*dr, dr)
num_increments = len(edges)-1
g = zeros([num_interior_particles, num_increments])
numberDensity = len(x)/S**2

# Compute pairwise correlation for each interior particle
for p in range(num_interior_particles):
index = interior_indices[p]
d = sqrt((x[index]-x)**2 + (y[index]-y)**2)
d[index] = 2*rMax

(result,bins) = histogram(d, bins=edges, normed=False, new=True)
g[p,:] = result/numberDensity

# Average g(r) for all interior particles and compute radii
g_average = zeros(num_increments)
for i in range(num_increments):
radii[i] = (edges[i] + edges[i+1])/2.
rOuter = edges[i+1]
rInner = edges[i]
g_average[i] = average(g[:,i])/(pi*(rOuter**2 - rInner**2))

return (g_average, radii, x[interior_indices], y[interior_indices])

#my code starts here
nlat=20	#HOW MANY DISKS ON THE SIDES OF THE BOX

r = 0.5	#DISK RADIUS

d = 2.0*r	#DIAMETER

s = 1.22	#GAPS, ALTER TO CHANGE DENSITY
S = s - 1

L = (nlat*d)+(nlat-1)*S	#LENGTH OF SIDES OF BOX
rMax = 5.0
dr=0.08 # Radius for increasing the size of the bins
S = L
area = L**2
pi = 3.141592654
diskarea= pi*(r**2) * nlat**2
density = diskarea/area
cpd = 0.9069
cpdf = (density / cpd)*100

print 'density = %f' % (density)
print 'which is %f percent of the closest packed density' % (cpdf)

scene.autoscale=0  #KEEP CAMERA STILL

scene = display(title='Monte Carlo Simulation of a Hard Disk Fluid',width=600, height=600,center=(L/2,L/2,0), background=(1,1,1))

#################################DRAW THE BOX#########################################

thk=0.3

s2 = L - thk

s3 = L + thk

wallR = box (pos=vector(L, L/2, 0), length=thk, height=s2, width=d, color = color.yellow)

wallL = box (pos=vector(0, L/2, 0), length=thk, height=s2, width=d, color = color.yellow)

wallB = box (pos=vector(L/2, 0, 0), length=s3, height=thk, width=d, color = color.yellow)

wallT = box (pos=vector(L/2, L, 0), length=s3, height=thk, width=d, color = color.yellow)

###################################################################################

#######CREATE SQUARE ARRAY OF PARTICLES (NLAT BY NLAT)############

alldisks = []

for x in range(0,nlat):

for y in range(0,nlat):

disk=sphere(pos=(x*s,y*s), radius = r, color=color.red)

alldisks.append(disk)

disk.id = (x,y)

####################################################

nstep=0.0				#SET FIRST STEP = 0

reject=0.0		#COUNT THE NUMBER OF TIMES A STEP IS REJECTED

maxpos=L		#DEFINE BOUNDARY CONDITIONS

while (nstep<5000):####HOW MANY STEPS?####

nstep=nstep+1

p = choice(alldisks) #CHOOSE A PARTICLE,P, AT RANDOM FROM THE ARRAY

pnewx=p.x + 0.9*(random() -0.5) #GENERATE A STEP IN X CHANGE COEFFICIENT OF RANDOM TO GET 30% ACCEPTANCE

pnewy=p.y + 0.9*(random() -0.5) #GENERATE A STEP IN Y

######PERIODIC BOUNDARY CONDITIONS############

if pnewx > L:			#IF X GOES BEYOND THE RIGHT WALL

pnewx=pnewx - L     #MOVE X POSITION TO OTHER SIDE

if pnewx < 0:		#IF X FALLS BELOW THE LEFT WALL

pnewx=pnewx + L #MOVE X POSITION TO OTHER SIDE

if pnewy > L:			#IF Y IS ABOVE BOX

pnewy=pnewy - L	    #MOVE TO THE BOTTOM

if pnewy < 0:		#IF Y BELOW

pnewy=pnewy + L #TO TOP

#################################################

rejected = False

for disk in alldisks:

distance=mag((vector(pnewx,pnewy))-disk.pos)#WORK OUT THE DISTANCE BETWEEN THE NEW POSITION FOR THE DISK AND ANY OTHER DISK

if (distance <= 2.0 * r) and p.id != disk.id: #if p overlaps any other disk

rejected = True

reject=reject+1 # COUNT HOW MANY TIMES THE STEP IS REJECTED

break

if not rejected:

p.x=pnewx	#

p.y=pnewy   #ACCEPT THE STEPS AS THE NEW POSITION FOR THE PARTICLE, P

if (nstep % 1800) == 0: #DUMP AFTER SO MANY STEPS
x_list = []
y_list = []
for disk in alldisks:
x_list.append(disk.x)
y_list.append(disk.y)
output = open('output''%s' % (nstep), 'w') #open file for outputting
(g, radii, newx, newy) = PairCorrelationFunction_2D(array(x_list), array(y_list), S, rMax, dr) #CALCULATE PAIR CORRELATION FUNCTION FOR EACH DUMP

for index in range(len(g)):
output.write( "%s       %s\n" % (radii[index], g[index]))

output.close()							#close file
continue

rejection = (reject/nstep) * 100 #percentage of steps rejected
acceptance = 100 - rejection
print 'acceptance rate was %f per cent' % (acceptance)

end = time.time()
elapsed= end - start
print "The simulation took", elapsed, "seconds to run"``````

## All 9 Replies

It seems, that most of the time is spent with this line:
distance=mag((vector(pnewx,pnewy))-disk.pos)

This line is repeated nlat*5000 times.

If you compute the vector before the for disk in alldisks loop, you gain some speed. In my measurement appr. 30%.

Some random remarks:
Do not use star import (from anything import *). You cannot know which object from which module comes. You could accidently shadow some objects too.
Do not import anything inside a function.
If you want to use psyco, don't use mainlevel code. Everything should be put into functions and classes.

Some words about the process profiling code.

To get profile data, you can use:
python -m profile script.py

If that does not give useful results (which is your case), you should but suspicious parts of the code into functions. In your case there are two points.:

The place where you set the p.x and p.y, because that is the point where the gui is updated, I think.
The other is the mentioned magnitude computation.

If i put the line outside the loop it only calculates the vector once, between one disk and another. It needs to loop through all the disks to ensure that none overlap.

Maybe I was not clear enough.

``````for disk in alldisks:
distance=mag((vector(pnewx,pnewy))-disk.pos)``````

use this

``````v=vector(pnewx,pnewy)
for disk in alldisks:
distance=mag(v-disk.pos)``````
commented: Very useful post to speed my program up! +1

ah of course!
Also I'm happy to run the program without the graphical parts,
But how would I create an object(disk) with the attributes in the above program myself?
I'm still a beginner.

If you ask, how to get rid of the visual module...
Well, I don't know.

Maybe numpy has some modules, which are modelling 2D objects.

Maybe the visual module can be told, not to display visuals:)
Try: scene.hide()

That's great! Shaved another 10 seconds off my benchmark time.
Many thanks

Some other improvements:

``````while (nstep<5000):####HOW MANY STEPS?####
nstep+=1
p = choice(alldisks) #CHOOSE A PARTICLE,P, AT RANDOM FROM THE ARRAY
pnewx=p.x + 0.9*(random()-0.5) #GENERATE A STEP IN X CHANGE COEFFICIENT OF RANDOM TO GET 30% ACCEPTANCE
pnewy=p.y + 0.9*(random()-0.5) #GENERATE A STEP IN Y

######PERIODIC BOUNDARY CONDITIONS############
if pnewx > L:		#IF X GOES BEYOND THE RIGHT WALL
pnewx-= L         #MOVE X POSITION TO OTHER SIDE
if pnewx < 0:		#IF X FALLS BELOW THE LEFT WALL
pnewx += L        #MOVE X POSITION TO OTHER SIDE
if pnewy > L:		#IF Y IS ABOVE BOX
pnewy-= L	      #MOVE TO THE BOTTOM
if pnewy < 0:		#IF Y BELOW
pnewy+=L          #TO TOP
#################################################

rejected = False
v=vector(pnewx,pnewy)
for disk in alldisks:
distance=mag(v-disk.pos)#WORK OUT THE DISTANCE BETWEEN THE NEW POSITION FOR THE DISK AND ANY OTHER DISK
if (distance <= 2.0 * r) and p is not disk: #if p overlaps any other disk
rejected = True
reject+=1 # COUNT HOW MANY TIMES THE STEP IS REJECTED
break

if not rejected:
p.x=pnewx	#
p.y=pnewy   #ACCEPT THE STEPS AS THE NEW POSITION FOR THE PARTICLE, P

if (nstep % 1800) == 0: #DUMP AFTER SO MANY STEPS
x_list = [disk.x for disk in alldisks]
y_list = [disk.y for disk in alldisks]
output = open('output''%s' % (nstep), 'w') #open file for outputting
(g, radii, newx, newy) = PairCorrelationFunction_2D(array(x_list), array(y_list), S, rMax, dr) #CALCULATE PAIR CORRELATION FUNCTION FOR EACH DUMP
lg=len(g)
for index in xrange(lg):
output.write( "%s       %s\n" % (radii[index], g[index]))
output.close()							#close file``````

Thanks, implementing most of that sped things up a little apart from the lines:

``````x_list = [disk.x for disk in alldisks]
y_list = [disk.y for disk in alldisks]``````

Changing this seemed to slow it down a couple of seconds.

Here is my full code though now:

``````from visual import *	#

from math import sqrt	#IMPORT LIBRARIES

from random import random, choice	#
from numpy import zeros, sqrt, where, pi, average, arange, histogram, array

import time

start = time.time()

def PairCorrelationFunction_2D(x,y,S,rMax,dr):
"""Compute the two-dimensional pair correlation function, also known
as the radial distribution function, for a set of circular particles
contained in a square region of a plane.  This simple function finds
reference particles such that a circle of radius rMax drawn around the
particle will fit entirely within the square, eliminating the need to
compensate for edge effects.  If no such particles exist, an error is
returned. Try a smaller rMax...or write some code to handle edge effects! ;)

Arguments:
x               an array of x positions of centers of particles
y               an array of y positions of centers of particles
S               length of each side of the square region of the plane
rMax            outer diameter of largest annulus
dr              increment for increasing radius of annulus

Returns a tuple: (g, radii, interior_x, interior_y)
g(r)            a numpy array containing the correlation function g(r)
radii           a numpy array containing the radii of the
annuli used to compute g(r)
interior_x      x coordinates of reference particles
interior_y      y coordinates of reference particles
"""
from numpy import zeros, sqrt, where, pi, average, arange, histogram
# Number of particles in ring/area of ring/number of reference particles/number density
# area of ring = pi*(r_outer**2 - r_inner**2)

# Find particles which are close enough to the box center that a circle of radius
# rMax will not cross any edge of the box

bools1 = x>1.1*rMax
bools2 = x<(S-1.1*rMax)
bools3 = y>rMax*1.1
bools4 = y<(S-rMax*1.1)
interior_indices, = where(bools1*bools2*bools3*bools4)
num_interior_particles = len(interior_indices)

if num_interior_particles < 1:
raise  RuntimeError ("No particles found for which a circle of radius rMax\
will lie entirely within a square of side length S.  Decrease rMax\
or increase the size of the square.")

edges = arange(0., rMax+1.1*dr, dr)
num_increments = len(edges)-1
g = zeros([num_interior_particles, num_increments])
numberDensity = len(x)/S**2

# Compute pairwise correlation for each interior particle
for p in range(num_interior_particles):
index = interior_indices[p]
d = sqrt((x[index]-x)**2 + (y[index]-y)**2)
d[index] = 2*rMax

(result,bins) = histogram(d, bins=edges, normed=False, new=True)
g[p,:] = result/numberDensity

# Average g(r) for all interior particles and compute radii
g_average = zeros(num_increments)
for i in range(num_increments):
radii[i] = (edges[i] + edges[i+1])/2.
rOuter = edges[i+1]
rInner = edges[i]
g_average[i] = average(g[:,i])/(pi*(rOuter**2 - rInner**2))

return (g_average, radii, x[interior_indices], y[interior_indices])

#my code starts here
nlat=20	#HOW MANY DISKS ON THE SIDES OF THE BOX

r = 0.5	#DISK RADIUS

d = 2.0*r	#DIAMETER

s = 1.22	#GAPS, ALTER TO CHANGE DENSITY
S = s - 1

L = (nlat*d)+(nlat-1)*S	#LENGTH OF SIDES OF BOX
rMax = 5
dr=0.1 # Radius for increasing the size of the bins
S = L
area = L**2
pi = 3.141592654
diskarea= pi*(r**2) * nlat**2
density = diskarea/area
cpd = 0.9069
cpdf = (density / cpd)*100

print 'density = %f' % (density)
print 'which is %f percent of the closest packed density' % (cpdf)

scene.hide()

#######CREATE SQUARE ARRAY OF PARTICLES (NLAT BY NLAT)############

alldisks = []

for x in range(0,nlat):

for y in range(0,nlat):

disk=sphere(pos=(x*s,y*s), radius = r, color=color.red)

alldisks.append(disk)

disk.id = (x,y)

####################################################

nstep=0.0				#SET FIRST STEP = 0

reject=0.0		#COUNT THE NUMBER OF TIMES A STEP IS REJECTED

maxpos=L		#DEFINE BOUNDARY CONDITIONS

print 'Running Simulation.....'

while (nstep<5000):####HOW MANY STEPS?####

nstep+=1

p = choice(alldisks) #CHOOSE A PARTICLE,P, AT RANDOM FROM THE ARRAY

pnewx=p.x + 0.9*(random() -0.5) #GENERATE A STEP IN X CHANGE COEFFICIENT OF RANDOM TO GET 30% ACCEPTANCE

pnewy=p.y + 0.9*(random() -0.5) #GENERATE A STEP IN Y

######PERIODIC BOUNDARY CONDITIONS############
if pnewx > L:		#IF X GOES BEYOND THE RIGHT WALL
pnewx-= L         #MOVE X POSITION TO OTHER SIDE
if pnewx < 0:		#IF X FALLS BELOW THE LEFT WALL
pnewx += L        #MOVE X POSITION TO OTHER SIDE
if pnewy > L:		#IF Y IS ABOVE BOX
pnewy-= L	      #MOVE TO THE BOTTOM
if pnewy < 0:		#IF Y BELOW
pnewy+=L          #TO TOP
#################################################

rejected = False
v=vector(pnewx,pnewy)

for disk in alldisks:

distance=mag(v-disk.pos)#WORK OUT THE DISTANCE BETWEEN THE NEW POSITION FOR THE DISK AND ANY OTHER DISK

if (distance <= 2.0 * r) and p.id != disk.id: #if p overlaps any other disk

rejected = True

reject+=1 # COUNT HOW MANY TIMES THE STEP IS REJECTED

break

if not rejected:

p.x=pnewx	#

p.y=pnewy   #ACCEPT THE STEPS AS THE NEW POSITION FOR THE PARTICLE, P

if (nstep % 1800) == 0: #DUMP AFTER SO MANY STEPS
x_list = []
y_list = []
for disk in alldisks:
x_list.append(disk.x)
y_list.append(disk.y)
output = open('output''%s' % (nstep), 'w') #open file for outputting
(g, radii, newx, newy) = PairCorrelationFunction_2D(array(x_list), array(y_list), S, rMax, dr) #CALCULATE PAIR CORRELATION FUNCTION FOR EACH DUMP
lg=len(g)

for index in range(lg):
output.write( "%s       %s\n" % (radii[index], g[index]))

output.close()							#close file
continue

rejection = (reject/nstep) * 100 #percentage of steps rejected
acceptance = 100 - rejection
print 'acceptance rate was %f per cent' % (acceptance)

end = time.time()
elapsed= end - start
print "The simulation took", elapsed, "seconds to run"``````
Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, learning, and sharing knowledge.