Ad:

Python Discussion Thread
Unsolved
Views: 669

Mar 10th, 2009

Python speed issue.

Expand Post »

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.

python Syntax (Toggle Plain Text)
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"
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"

Last edited by breatheasier; Mar 10th, 2009 at 6:17 pm. Reason: removed some pointless whitespace

Similar Threads

Speeding up Python in Python
Dictionary Sorting? in Python
Non-programmer tutorial, stuck at list exercise in Python

breatheasier

Reputation Points: 11

Solved Threads: 2

Light Poster

Offline

46 posts
since Feb 2009

Permalink

Mar 10th, 2009

Re: Python speed issue.

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

Reputation Points: 56

Solved Threads: 65

Posting Whiz in Training

Offline

242 posts
since Jun 2008

Permalink

Mar 10th, 2009

Re: Python speed issue.

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.

slate

Reputation Points: 56

Solved Threads: 65

Posting Whiz in Training

Offline

242 posts
since Jun 2008

Permalink

Mar 10th, 2009

Re: Python speed issue.

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.

breatheasier

Reputation Points: 11

Solved Threads: 2

Light Poster

Offline

46 posts
since Feb 2009

Permalink

Mar 10th, 2009

Re: Python speed issue.

Maybe I was not clear enough.

instead of this:

Python Syntax (Toggle Plain Text)
    for disk in alldisks:
        distance=mag((vector(pnewx,pnewy))-disk.pos)
    for disk in alldisks:
        distance=mag((vector(pnewx,pnewy))-disk.pos)

use this

Python Syntax (Toggle Plain Text)
   v=vector(pnewx,pnewy)
   for disk in alldisks:
        distance=mag(v-disk.pos)
   v=vector(pnewx,pnewy)
   for disk in alldisks:
        distance=mag(v-disk.pos)

slate

Reputation Points: 56

Solved Threads: 65

Posting Whiz in Training

Offline

242 posts
since Jun 2008

Permalink

Mar 10th, 2009

Re: Python speed issue.

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.

breatheasier

Reputation Points: 11

Solved Threads: 2

Light Poster

Offline

46 posts
since Feb 2009

Permalink

Mar 10th, 2009

Re: Python speed issue.

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

slate

Reputation Points: 56

Solved Threads: 65

Posting Whiz in Training

Offline

242 posts
since Jun 2008

Permalink

Mar 10th, 2009

Re: Python speed issue.

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

breatheasier

Reputation Points: 11

Solved Threads: 2

Light Poster

Offline

46 posts
since Feb 2009

Permalink

Mar 10th, 2009

Re: Python speed issue.

Some other improvements:

Python Syntax (Toggle Plain Text)
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
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

slate

Reputation Points: 56

Solved Threads: 65

Posting Whiz in Training

Offline

242 posts
since Jun 2008

Permalink

Mar 10th, 2009

Re: Python speed issue.

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

Python Syntax (Toggle Plain Text)
        x_list = [disk.x for disk in alldisks]
        y_list = [disk.y for disk in alldisks]
        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:

python Syntax (Toggle Plain Text)
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"
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"

Last edited by breatheasier; Mar 10th, 2009 at 9:38 pm.

breatheasier

Reputation Points: 11

Solved Threads: 2

Light Poster

Offline

46 posts
since Feb 2009

This thread is more than three months old

No one has posted to this discussion for at least three months. Please let old threads die and do not reply to them unless you feel you have something new and valuable to contribute that absolutely must be added to make the discussion complete. Otherwise, please start a new thread in this forum instead.

Previous Thread in Python Forum Timeline: Error?

Next Thread in Python Forum Timeline: Question

DaniWeb Message

Cancel Changes

User Name
Password