Counting Neighbors (Why scipy.signal.convolve2D so fast?)

Question

Here is my Python implementation of counting neighbours of Game of Life with radius as parameter.

def neighbors_count(n2d_array, radii=1):
    assert n2d_array.ndim == 2
    row_len, col_len = n2d_array.shape
    nbrs_count = np.zeros_like(n2d_array)
    for row_idx, row_val in enumerate(n2d_array):
         for col_idx, col_val in enumerate(row_val):
             start_row = 0 if (row_idx-radii) < 0 else (row_idx-radii)
             end_row = row_len if (row_idx+radii+1) > row_len else (row_idx+radii+1)
             start_col = 0 if (col_idx-radii) < 0 else (col_idx-radii)
             end_col = row_len if (col_idx+radii+1) > row_len else (col_idx+radii+1)
             neighbor = 0
             for block_row_idx in np.arange(start_row, end_row):
                 for block_col_idx in np.arange(start_col, end_col):
                     neighbor += n2d_array[block_row_idx, block_col_idx]
             nbrs_count[row_idx, col_idx] = neighbor - n2d_array[row_idx, col_idx]
     return nbrs_count

I found out that my implementation is very slow compared to scipy.signal.convolve2d:

def neighbors_count2(n2d_array, radii=1):
    from scipy.signal import convolve2d
    diameter = 2 * radii + 1
    n2d_array = n2d_array.astype(bool)
    nbrs_count = convolve2d(n2d_array, np.ones((diameter, diameter)),
                            mode='same', boundary='fill') - n2d_array
    return nbrs_count

Here is %timeit result in my computer:

%timeit -n 10 neighbors_count(np.random.randint(2, size=(100,100)))
10 loops, best of 3: 232 ms per loop

%timeit -n 10 neighbors_count2(np.random.randint(2, size=(100,100)))
10 loops, best of 3: 963 µs per loop

How to improve/vectorize my code so it can run faster than scipy.signal.convolve2d?

Answer 1

You can change your algorithmic approach to improve your speed.

Currently:

• look at every cell
• lookup every neighbor-cell
• write the number of neighbors

My proposition: You start with a zero-ed nbrs_count array and look at every cell. If it's occupied, increase the nbrs_count of all neighboring cells (you will get a huge speedup if the array is mostly empty).

To prevent all your conditional statements, you can simply use a try: ... except: block, as suggested by @JoeWallis

---

Here is an implementation using my propositions:

import numpy as np

def neighbors_count(n2d_array, radii=1):
    assert n2d_array.ndim == 2

    nbrs_count = np.zeros_like(n2d_array)

    # array of adjacents cells
    adjacents = []
    for i in range(-radii, radii + 1):
        for j in range(-radii, radii + 1):
            if j != 0 or i != 0:
                adjacents.append((i,j))

    for (row_idx, col_idx), value in np.ndenumerate(n2d_array):
        if value:
            for i, j in adjacents:
                try:
                    if row_idx + i >= 0 and col_idx + j >= 0:
                        # because a negative index doesn't fail
                        nbrs_count[row_idx + i, col_idx + j] += 1
                except IndexError:
                    # We are outside the array
                    pass

    return nbrs_count

This solution is about 5 times faster than the original code (which is still way slower than scipy)

Answer 2

So to get the questions out of the way first.

Why scipy.signal.convolve2D so fast?

SciPy and Co. all are programmed in C, with a Python interface. As almost everyone knows interpreted languages are slow compared to compiled languages for the most part and will explain the difference in speed.
For a very small seemingly unbiased comparison of C vs Python, you can look at the Julia home page. It says Python can range from 15 to 30 times slower than C, again over a small set of functions and excluding anomaly's.

How to improve/vectorize my code so it can run faster than scipy.signal.convolve2d?

Unless you want a complete re-write of your entire game of life in say C/C++, good luck.

---

• It's best to not use assert, as they get ignored when you call the script with the -O flag.

• Limit the length of your lines. 79 is the common maximum character width in Python.

• Python prefers to have white-space over clumped up code, so when using operators try to leave a space either side of them. 1-1 looks like a variable not like 1 - 1. This is worse since 1 can look like l in some typefaces.

• You should change the turnery operators:

  • Expand them over multiple lines, to increase readability;

  • Change them to use a function, so you don't calculate row_idx - radii twice; or,

  • Use min and max. Which gains the benefits of both the above.

• You assume that you have \$O(1)\$ writes, neighbour has \$O(3^2\times\text{n2d_array})\$ reads and writes. You are computing and writing to start_row and end_row way more then you need to, try putting them in the above for loop.

• Your sum for neighbour can be written nicer using a comprehension and sum, but you may loose some performance.

• Some of your variable names should probably change, idx can just be index, and end_row seems more like a function name than the row_end.

I would write your code as the following:

def neighbours_count(array, radii=1):
    if array.ndim != 2:
        raise Exception # Change this to a better Exception

    row_size, column_size = array.shape
    neighbours_amount = np.zeros_like(array)
    for row_index, row_value in enumerate(array):
        row_start = max(row_index - radii, 0)
        row_end = min(row_index + radii + 1, row_size)
        for column_index, column_value in enumerate(row_value):
            column_start = max(column_index - radii, 0)
            column_end = min(column_index + radii + 1, column_size)

            neighbours_amount[row_index, column_index] = sum(
                array[i, j]
                for i in np.arange(start_row, end_row)
                for j in np.arange(start_column, end_column)
            ) - array[row_index, column_index]

    return neighbours_amount

As there is no feasible way to make Python outperform scipy.signal.convolve2d I didn't improve the performance.

Instead I would recommend either:

1. Don't re-invent the wheel and use scipy.signal.convolve2d, or;
2. Use @oliverpool's answer.