# 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. – twohundredping Jun 10 '15 at 22:43
• Ok, I see that it can be separated into two vectors. – inyoot Jun 10 '15 at 23:44
• How is it implemented in C? The same loop structure as I wrote in Python? – inyoot Oct 24 '15 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. – Peilonrayz Oct 24 '15 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. – inyoot Oct 26 '15 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 – oliverpool Oct 26 '15 at 10:26
• I've implemented my proposition: please have a look at my updated answer – oliverpool Oct 26 '15 at 11:50
• Well, trying to beat C using Python seems quite ambitious... – oliverpool Oct 26 '15 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. – Peilonrayz Oct 30 '15 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/… – holroy Oct 27 '15 at 23:38