Have something to say? Contribute New Article Reply to this Article
Python speed issue.
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])
radii = zeros(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"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.
instead of this:
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)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()
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 fileThanks, 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])
radii = zeros(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"This article has been dead for over three months
You
© 2012 DaniWeb® LLC
