Python speed issue.

[breatheasier]

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"

[slate]

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.

[slate]

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.

[breatheasier]

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.

[slate]

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)

[breatheasier]

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.

[slate]

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()

[breatheasier]

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

[slate]

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

[breatheasier]

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])
    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"