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

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):
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?

• See separable filters for now. I am sure someone with more experience in Python will help you with the rest. Commented Jun 10, 2015 at 22:43
• Ok, I see that it can be separated into two vectors. Commented Jun 10, 2015 at 23:44
• How is it implemented in C? The same loop structure as I wrote in Python? Commented Oct 24, 2015 at 19:51
• I deleted my comment as I found the Python interface, but it's written in C. You can check the source... which links to a .py file that then links to the C file. Commented Oct 24, 2015 at 19:54

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

assert n2d_array.ndim == 2

nbrs_count = np.zeros_like(n2d_array)

if j != 0 or i != 0:

for (row_idx, col_idx), value in np.ndenumerate(n2d_array):
if value:
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)

• I don't see that will increase speed because it will read once and write many times instead of read many and write once. Writing is slower than reading. Commented Oct 26, 2015 at 4:23
• You are already writing many times your neighbor variable, I don't think that writing nbrs_count[row_idx, col_idx] will be so much slower Commented Oct 26, 2015 at 10:26
• I've implemented my proposition: please have a look at my updated answer Commented Oct 26, 2015 at 11:50
• Well, trying to beat C using Python seems quite ambitious... Commented Oct 26, 2015 at 13:20
• I was thinking about this, rather than extending the array, you could do a try except. nbrs_count = np.zeros_like(array), and try: nbrs_count[row_idx + i, col_idx + j] += 1. This would just make the entire function easier to read, as row_idx + i + radii is weird without the comment. And nbrs_count = np.delete(nbrs...) will not work with radii > 1. Commented Oct 30, 2015 at 12:30

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.

1. Don't re-invent the wheel and use scipy.signal.convolve2d, or;
• assert is ignored when -O is used in Python 2 or Python 3, according to docs.python.org/2/reference/… and docs.python.org/3/reference/… Commented Oct 27, 2015 at 23:38