Lattice lotka-volterra 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 np
import matplotlib.pyplot as plt
import random
import time
import 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 matrix

def mort(matrix,d):
    RN = np.random.random((n,n))
    casualty = (RN &lt; d)
    matrix[...][casualty] = 0 
    return matrix

def 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&lt;0 or i+p&gt;n-1) or (j+q&lt;0 or j+q&gt;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] &gt; np.random.random():
                matrix[i+p,j+q] = matrix[i,j]
    return matrix

def 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 rate
spp = 10
B = [0.57710226,0.57722738,0.577233935,0.577120507,0.576885731,0.576528293,
     0.576046935,0.575440457,0.574707715,0.573847625]
dr = 0.3

#initial density
initD = 0.005

#mortality rate
m = 0.1

#matrix size
n = 1000

#lattice system is initially empty
latsys = np.zeros((n,n))

#randomly place the species in the matrix
for i in range(1,spp+1):
    dispose(latsys,initD,spp,i)

#reproduction and mortality
for 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-t0
print 'Total simulation run time =', total4, 'second(s).'

#for graphing snapshots
plt.imshow(latsys, cmap=plt.cm.spectral, interpolation='nearest')
plt.colorbar()
plt.show()
There's quite a lot of code here, so for now I'm just going to look at the dispose function. You'll see that there's more than enough to say about this function for one answer.

1. Review


There's no docstring. What does this function do? How do I call it? What values do I pass for parameters? What does it return? It's very hard to maintain code when you don't know what it is supposed to do. This doesn't have to be long-winded or laborious: a quick explanation like the following is perfectly adequate.

"""Choose n*n*dens positions in matrix with the value 0, uniformly at
random, and set them to sp. Return matrix."""

There's no good reason to abbreviate the parameter to dens. It only appears once in the function body, and we're not running out of letters, so call it density.
The function takes a parameter s but this is not used. Better to omit it.
The function uses the global variable n. It would be better to pass this in as a parameter, or even better, deduce it from the size of matrix. Presumably the intention is that matrix is always an n × n array, in which case you should write matrix.size instead of n * n.
The name dispose is vague and does not convey much to me. I think the name populate would be better.
The function always returns matrix, which is useless, as the caller already knows about matrix.
The name matrix is slightly misleading, as Numpy has a matrix class but this parameter does not belong to it! It's actually an ndarray. I would call the parameter array or maybe just a.
The name sp is rather too specific. Presumably it is short for "species", but really nothing about this function depends on this value representing a species. You can imagine other types of lattice simulation where the array entries represent molecules, grains of sand, people, etc. So for generality I'd name this parameter value.
The value 0 is arbitrary and so for extra generality it should be a function parameter (taking the default value 0). Imagine having a "speciation" operation in which you replace some fraction of a species with a descendant. It would be nice to be able to reuse this code.
If dens is greater than 1, then limit will be greater than the number of elements in the array, and so the function will run forever.

Even if dens is less than 1, there might still not be enough positions in the array with value 0, and so the function will run forever.

Perhaps you plan not to fill very much of the array, but even so it is a good idea to write robust code in case you change your plans or want to reuse this code in other programs.
The way that the positions are found is to repeatedly generate random coordinates and check to see if the position is 0. This is inefficient. Imagine that dens is 1, so that all positions must be set to sp. How long will this take, on average? This is the well known "coupon collector's problem" so the expected number of coordinates to be randomly generated is \$ n^2 \log n^2 \$.

Perhaps you plan not to fill very much of the array, but again, it is a good idea to write code with the right asymptotic complexity, in case you change your plans or want to reuse this code in other programs.

So it would be better to find the positions containing 0 (using numpy.flatnonzero), then choose a random sample of those positions without replacement (using numpy.random.choice), and then set all the positions.


2. Performance

Putting all this together, I'd write the function like this:

def populate(array, density, value, original=0):
    """Choose array.size * density positions in array that are equal to
    original, uniformly at random, and set them to value.

    """
    a = array.reshape(-1)
    count = int(a.size * density)
    a[np.random.choice(np.flatnonzero(a == original), count, False)] = value


Note the use numpy.reshape to create a flattened view of array, so that we can use numpy.random.choice, which requires a 1-D array. We then use this view to update the original array.

However, this has unacceptably poor performance compared to dispose when density is small:

&gt;&gt;&gt; a = np.zeros((1000, 1000))
&gt;&gt;&gt; timeit(lambda:dispose(a, 0.01, None, 1), number=1)
0.0184222329990007
&gt;&gt;&gt; timeit(lambda:populate(a, 0.01, 1), number=1)
0.2894051409966778


The reason for this poor performance is that numpy.random.choice computes a sample without replacement by permuting the whole array and using an initial prefix of the permutation. This is expensively wasteful when only a small part of the permutation is needed.

So what about sampling with replacement instead (which is fast), eliminating duplicates using numpy.unique and looping until we have enough choices (as in the original post)? Here's the implementation:

def populate(array, density, value, original=0):
    """Choose array.size * density positions in array that are equal to
    original, uniformly at random, and set them to value.

    """
    a = array.reshape(-1)
    count = int(a.size * density)
    while count:
        p = np.random.choice(a.size, count)
        p = np.unique(p[a[p] == original])
        a[p] = value
        count -= p.size


This is adequately fast: it's about ten times as fast as the original code, and more than a hundred times as fast as my first attempt above:

&gt;&gt;&gt; a = np.zeros((1000,1000))
&gt;&gt;&gt; timeit(lambda:populate(a, 0.01, 1), number=1)
0.0019981740042567253


However, it's not robust: it suffers from the same problems as the original code, discussed in §1.10 and §1.11 above. What we really need is a hybrid implementation: that tries the approach that's fast for small densities, but measures how well it's doing, and switches to the slow but sure approach for large densities. For example:

def populate(array, density, value, original=0):
    """Choose array.size * density positions in array that are equal to
    original, uniformly at random, and set them to value.

    """
    a = array.reshape(-1)
    count = int(a.size * density)
    while count:
        p = np.random.choice(a.size, count)
        p = np.unique(p[a[p] == original])
        a[p] = value
        count -= p.size
        if count &gt; p.size * 10:
            # Not making much progress: switch to robust but slow algorithm
            a[np.random.choice(np.flatnonzero(a == original), count, False)] = value
            return


This remains very fast when density is small:

&gt;&gt;&gt; a = np.zeros((1000, 1000))
&gt;&gt;&gt; timeit(lambda:populate(a, 0.01, 1), number=1)
0.0022136839979793876


but still completes in a reasonable time when density is large:

&gt;&gt;&gt; a = np.zeros((1000, 1000))
&gt;&gt;&gt; timeit(lambda:populate(a, 0.99, 1), number=1)
0.3501391080208123
&gt;&gt;&gt; a = np.zeros((1000, 1000))
&gt;&gt;&gt; timeit(lambda:dispose(a, 0.99, None, 1), number=1)
5.953548110992415


and doesn't go into an infinite loop if there aren't enough spaces remaining:

&gt;&gt;&gt; timeit(lambda:populate(a, 1.1, 1), number=1)
Traceback (most recent call last):
  File "&lt;stdin&gt;", line 1, in &lt;module&gt;
  ...
  File "mtrand.pyx", line 1106, in mtrand.RandomState.choice (numpy/random/mtrand/mtrand.c:11223)
ValueError: Cannot take a larger sample than population when 'replace=False'