Can anyone comment on the efficiency of my code? I am new with Python using NumPy/SciPy and I am trying to program a lattice biological model where a certain (one and only one) species (represented by a nonzero element) occupies a lattice square (represented in the code by a 2D array). This species can reproduce to an adjacent and empty lattice square (represented by a zero element). (Quite similar to a lattice gas model.)The part of the code that takes the longest time is the reproduction part (def repro) because it has to check whether a site is already chosen every time a random site is generated to avoid repetition (the first while loop in def repro). Can anyone suggest a faster way to do this?I am thinking of creating a collection of indices (2D index) where I can choose a random element and delete it at the same time so that the while loop won't be necessary. How can I do this?import numpy as npimport matplotlib.pyplot as pltimport randomimport timeimport csv open('p14.4.csv','w').close()t0 = time.time()def dispose(matrix,dens,s,sp):    limit = int(n*n*dens)    for num in range(limit):        i = np.random.randint(n)        j = np.random.randint(n)                while matrix[i,j] != 0:            i = np.random.randint(n)            j = np.random.randint(n)        matrix[i,j] = sp    return matrixdef mort(matrix,d):    RN = np.random.random((n,n))    casualty = (RN < d)    matrix[...][casualty] = 0     return matrixdef reprod(matrix,b):    p = np.random.randint(-1,2)    q = np.random.randint(-1,2)    idx = np.random.randint(-1,0,(n,n))    jdx = np.random.randint(-1,0,(n,n))    for k in range(n*n):        i = np.random.randint(n)        j = np.random.randint(n)        while (idx[i,j] == i and jdx[i,j] == j):            i = np.random.randint(n)            j = np.random.randint(n)        idx[i,j], jdx[i,j] = i, j        if matrix[i,j] != 0:            spp = int(matrix[i,j])            while ((i+p<0 or i+p>n-1) or (j+q<0 or j+q>n-1)):                p = np.random.randint(-1,2)                q = np.random.randint(-1,2)            if matrix[i+p,j+q] == 0 and b[spp-1] > np.random.random():                matrix[i+p,j+q] = matrix[i,j]    return matrixdef csv_writer(data,path,line):    with open(path, a) as csv_file:        writer = csv.writer(csv_file, delimiter=',')        writer.writerow(data)        csv_file.close()    return path#birth ratespp = 10B = [0.57710226,0.57722738,0.577233935,0.577120507,0.576885731,0.576528293,     0.576046935,0.575440457,0.574707715,0.573847625]dr = 0.3#initial densityinitD = 0.005#mortality ratem = 0.1#matrix sizen = 1000#lattice system is initially emptylatsys = np.zeros((n,n))#randomly place the species in the matrixfor i in range(1,spp+1):    dispose(latsys,initD,spp,i)#reproduction and mortalityfor t in range(1):    print t    x=latsys        #reproduction    reprod(latsys,B)    t3 = time.time()    total3 = t3-t0        print total3    #mortality    mort(latsys,dr)    S1 = float(np.count_nonzero(x==1))/(n*n)    S2 = float(np.count_nonzero(x==2))/(n*n)    S3 = float(np.count_nonzero(x==3))/(n*n)    S4 = float(np.count_nonzero(x==4))/(n*n)    S5 = float(np.count_nonzero(x==5))/(n*n)    S6 = float(np.count_nonzero(x==6))/(n*n)    S7 = float(np.count_nonzero(x==7))/(n*n)    S8 = float(np.count_nonzero(x==8))/(n*n)    S9 = float(np.count_nonzero(x==9))/(n*n)    S10 = float(np.count_nonzero(x==10))/(n*n)        print S1,S2,S3,S4,S5,S6,S7,S8,S9,S10    R = [t,S1,S2,S3,S4,S5,S6,S7,S8,S9,S10]    path = p14.4.csv    csv_writer(R,path,t)SPCS = [S1,S2,S3,S4,S5,S6,S7,S8,S9,S10]SURVIVORS = np.count_nonzero(SPCS)print The number of surviving species is , SURVIVORS  t4 = time.time()total4 = t4-t0print 'Total simulation run time =', total4, 'second(s).'#for graphing snapshotsplt.imshow(latsys, cmap=plt.cm.spectral, interpolation='nearest')plt.colorbar()plt.show()
Lattice lotka-volterra
python;numpy;scipy;simulation