Code speedup tips

[Carl Banks]

Sean Richards wrote:
> Hi,

Hello.


[snip]
> I can see that the nested loops are the bottleneck but what better
> alternatives does Python have for iterating over a 2d array.

The nested loops are the bottleneck; however, Numeric does have
alternatives.  In fact, Numeric was invented so that you don't have to
iterate over arrays (at the Python level, anyways).  Iterating over a
Numeric array is almost always the wrong way to do it.

What you want to do is take advantage of Numeric's built-in iteration,
which is much faster.  I'm sure you're aware that adding two arrays
adds all the elements in them.  You can apply the automaton in the
same way.

Consider the simple automaton rule here.  A cell's value in the
following iteration depends on the values of four pixels from the
current iteration: the cells immediately above, below, left of, and
right of the cell in question.  A simple way to apply this rule is to
add the four pixel values, and the resulting pixel is 1 if the sum is
nonzero.  We can add the pixel values 

Now, you might be thinking, "But Numeric operations only work with
corresponding entries.  I can add a[3,2] to b[3,2], but not a[3,2] to
a[3,4].  How can I do this quickly with Numeric?"

The answer is to use slicing.  What you want to do is take various
slices from the lattice, each of them of the same size, but offset so
that when you add them, the proper cells line up.  

To exemplify: suppose you have an array A of dimension 10x10.  Let B
and C be the following slices of A:

    B = A[:-2,:-2]
    C = A[2:,2:]

Now, let D = B + C.  Both B and C are 8x8 arrays, so this is possible.
Consider what the result D[3,3] will be.  It would be B[3,3]+C[3,3],
of course.  B[3,3] is A[3,3], while C[3,3] is A[5,5].  The end result
is that D[3,3] == A[3,3]+A[5,5].  Using this technique, the loops can
be gotten rid of.

Understand?  It's kind of hard to put your head around, I imagine,
especially if you don't understand slicing.

Here is how I rewrote CA_Function.  There are some details in it I've
left unexplained; feel free to ask for a clarification:


    def CA_Function(Length,Steps):
        # Create the two lattices
        current = zeros((Length,Length),Int)
        workspace = zeros((Length-2,Length-2),Int)

        # Create slices of the current lattice
        # We don't have to do this every step because slicing shares memory
        left = current[0:-2,1:-1]
        up = current[1:-1,0:-2]
        right = current[2:,1:-1]
        down = current[1:-1,2:]

        # Set the start cell to black
        current[Length/2,Length/2] = 1

        # Apply the rule
        for step in xrange(Steps):
            # In the workspace, add the shifted lattices
            # This uses a bunch of += because it's faster
            workspace[:] = left
            workspace[:] += up
            workspace[:] += right
            workspace[:] += down
            # Set the current lattice
            current[1:-1,1:-1] = where(workspace,1,0)


-- 
CARL BANKS