How to run my python code in Cython?

Question

I have currently made quite a large code in python and when I run it, it takes about 3 minutes for it to make the full calculation. Eventually, I want to increase my N to about 400 and change my m in the for loop to an even larger number - This would probably take hours to calculate which I want to cut down.

It's steps 1-6 that take a long time.

When attempting to run this with cython (I.E. importing pyximport then importing my file) I get the following error FDC.pyx:49:19: 'range' not a valid cython language construct and FDC.pyx:49:19: 'range' not a valid cython attribute or is being used incorrectly

from physics import *
from operator import add, sub
import pylab




################ PRODUCING CHARGES AT RANDOM IN r #############

N=11 #Number of point charges

x = zeros(N,float) #grid
y = zeros(N,float)

i=0
while i < N: #code to produce values of x and y within r 
    x[i] = random.uniform(0,1)
    y[i] = random.uniform(0,1)
    if x[i] ** 2 + y[i] ** 2 <= 1:
        i+=1


print x, y

def r(x,y): #distance between particles
    return sqrt(x**2 + y**2)   

o = 0; k = 0; W=0        #sum of energy for initial charges 
for o in range(0, N):
    for k in range(0, N):
        if o==k:
            continue
        xdist=x[o]-x[k]
        ydist=y[o]-y[k]
        W+= 0.5/(r(xdist,ydist))

print "Initial Energy:", W


##################### STEPS 1-6 ######################

d=0.01 #fixed change in length
charge=(x,y)
l=0; m=0; n=0
prevsW = 0.
T=100
for q in range(0,100):
    T=0.9*T
    for m in range(0, 4000): #steps 1 - 6 in notes looped over

        xRef = random.randint(0,1)      #Choosing x or y
        yRef = random.randint(0,N-1)      #choosing the element of xRef
        j = charge[xRef][yRef]           #Chooses specific axis of a charge and stores it as 'j'

        prevops = None #assigns prevops as having no variable
        while True: #code to randomly change charge positions and ensure they do not leave the disc
            ops =(add, sub); op=random.choice(ops)
            tempJ = op(j, d)
            #print xRef, yRef, n, tempJ
            charge[xRef][yRef] = tempJ
            ret = r(charge[0][yRef],charge[1][yRef])
            if ret<=1.0:
                j=tempJ
                #print "working", n
                break

            elif prevops != ops and prevops != None: #!= is 'not equal to' so that if both addition and subtraction operations dont work the code breaks 
                break

            prevops = ops #####

        o = 0; k = 0; sW=0        #New energy with altered x coordinate
        for o in range(0, N):
            for k in range(0, N):
                if o==k:
                    continue
                xdist = x[o] - x[k]
                ydist = y[o] - y[k]
                sW+=0.5/(r( xdist , ydist )) 


        difference = sW - prevsW
        prevsW = sW

        #Conditions:
        p=0

        if difference < 0: #accept change
            charge[xRef][yRef] = j
            #print 'step 5'
        randomnum = random.uniform(0,1) #r
        if difference > 0: #acceptance with a probability 
            p = exp( -difference / T )
            #print 'step 6', p
            if randomnum >= p:
                charge[xRef][yRef] = op(tempJ, -d) #revert coordinate to original if r>p
                #print charge[xRef][yRef], 'r>p'

        #print m, charge, difference   

o = 0; k = 0; DW=0        #sum of energy for initial charges 
for o in range(0, N):
    for k in range(0, N):
        if o==k:
            continue
        xdist=x[o]-x[k]
        ydist=y[o]-y[k]
        DW+= 0.5/(r(xdist,ydist))

print charge
print 'Final Energy:', DW

################### plotting circle ###################

# use radians instead of degrees
list_radians = [0]


for i in range(0,360):
    float_div = 180.0/(i+1)
    list_radians.append(pi/float_div)

# list of coordinates for each point
list_x2_axis = []
list_y2_axis = []

# calculate coordinates 
# and append to above list
for a in list_radians:
    list_x2_axis.append(cos(a))
    list_y2_axis.append(sin(a))



# plot the coordinates
pylab.plot(list_x2_axis,list_y2_axis,c='r')

########################################################
pylab.title('Distribution of Charges on a Disc') 
pylab.scatter(x,y)
    pylab.show()

Answer 1

What is taking time seems to be this:

for q in range(0,100):
    ...
    for m in range(0, 4000): #steps 1 - 6 in notes looped over
        while True: #code to randomly change charge positions and ensure they do not leave the disc
            ....
        for o in range(0, N):      # <----- N will be brought up to 400
            for k in range(0, N):
                ....
            ....            
        ....
    ....

100 x 4000 x (while loop) + 100 x 4000 x 400 x 400 = [400,000 x while loop] + [64,000,000,000]

Before looking into a faster language, maybe there is a better way to build your simulation?

Other than that, you will likely have immediate performance gains if you: - shift to numpy arrays i/o python lists. - Use xrange i/o range

[edit to try to answer question in the comments]:

import numpy as np, random
N=11 #Number of point charges

x = np.random.uniform(0,1,N)  
y = np.random.uniform(0,1,N)
z = np.zeros(N)
z = np.sqrt(x**2 + y**2)   # <--- this could maybe replace r(x,y) (called quite often in your code)
print x, y, z

You also could look into all the variables that are assigned or recalculated many many times inside your main loop (the one described above), and pull all that outside the loop(s) so it is not repeatedly assigned or recalculated.

for instance,

ops =(add, sub); op=random.choice(ops)

maybe could be replaced by

ops = random.choice(add, sub)

Lastly, and here I am out on a limb because I've never used it myself, but it might be a little bit simpler for you to use a package like Numba or Jit as opposed to cython; they allow you to decorate a critical part of your code and have it precompiled prior to execution, with none or very minor changes.