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

• See separable filters for now. I am sure someone with more experience in Python will help you with the rest. Jun 10, 2015 at 22:43
• Ok, I see that it can be separated into two vectors. Jun 10, 2015 at 23:44
• How is it implemented in C? The same loop structure as I wrote in Python? 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. Oct 24, 2015 at 19:54

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

assert n2d_array.ndim == 2

nbrs_count = np.zeros_like(n2d_array)

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

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)

• 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. 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 Oct 26, 2015 at 10:26
• I've implemented my proposition: please have a look at my updated answer Oct 26, 2015 at 11:50
• Well, trying to beat C using Python seems quite ambitious... 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. 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.

Instead I would recommend either:

1. Don't re-invent the wheel and use scipy.signal.convolve2d, or;
2. Use @oliverpool's answer.
• 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/… Oct 27, 2015 at 23:38
• How about using tree structure like hashlife? Will it increase performance? Oct 27, 2015 at 23:46
• @inyoot Python is a tool for fast development time. It is not meant to execute fast. Also you only provided one function, I cannot review or make implementations of code I have no clue about. Oct 27, 2015 at 23:52
• I just wondering is the speed-up mainly because of compiled language or using better algorithm and data structure? Oct 27, 2015 at 23:59
• @inyoot I don't know enough about the algorithm they use, but C can run 15-30 times faster than Python, also you can do some fancy stuff in it to make it faster. If you just want execution speed then compiled languages are generally faster than interpreted languages. I would assume if you wrote oliverpool's code in C it would have roughly the same speed. Oct 28, 2015 at 0:07