Speed of Python vs Julia for Game of Life. Why is Julia slow?

To compare the speed of Julia to Python+numba, I implemented Game of Life in both languages.

For anyone who is not familiar with Game of Life: Start with a matrix of boolean entries, encoding whether a cell is alive or dead. Then, at every iteration, update the state of a cell based upon its value and the eight neighboring values.

1. If the cell is alive, and exactly two or three neighboring cells are alive, the cell stays alive.
2. If the cell is alive and fewer than two or more than three neighboring cells are alive, the cell dies.
3. If a cell is dead, and exactly three neighbors are alive, it becomes alive.
4. If none of the above applies, the cell keeps it state.

I chose a special initial condition called Acorn, just to make sure the population stays alive sufficiently long.

Here is the reference Python+numba code with the respective timing and output.

import numpy as np
import matplotlib.pyplot as plt
from time import time
from numba import njit

def acorn(n): # Initial condition
state = np.zeros((n,n))
mid = n//2
state[mid,mid-2:mid+5] = 1
state[mid,mid] = 0
state[mid,mid+1] = 0
state[mid+1,mid+1] = 1
state[mid+2,mid-1] = 1
return state>0

@njit
def update(s):
dimx,dimy = s.shape
nexts = np.copy(s)
for i in range(1,dimx-1):
for j in range(1,dimy-1):
count = -1 if s[i,j] else 0
for offi in [-1,0,+1]:
for offj in [-1,0,+1]:
if s[i+offi,j+offj]:
count += 1

if s[i,j]:
if count < 2 or count > 3:
nexts[i,j] = False
else:
if count ==3:
nexts[i,j] = True
return nexts

def run(state,n):
for i in range(n):
state = update(state)
return state

n = 100
nruns = 1000
state = acorn(n)
t = time()
s = run(state,nruns)
dt = time() -t
plt.pcolormesh(s)
plt.grid("on")
print("Elapsed time {:.5f} seconds".format(dt))
print("One run in {:.5f} seconds".format(dt/n))
print("{} runs per seconds".format(int(n/dt)))
plt.savefig("out.png")


Output (after a precompilation run):

Elapsed time 0.07883 seconds
One run in 0.00079 seconds
1268 runs per seconds


This seems reasonably fast on my Intel i5 notebook with 12gb of RAM and Python 3.8 and recent numpy + numba versions.

Next is the Julia code. I basically translated the Python code to Julia using my limited Julia knowledge.

using PyPlot

function acorn(n)
mid = Int(n/2)+1
state = zeros(n,n)
state[mid,mid-2:mid+4] .= 1
state[mid,mid] = 0
state[mid,mid+1] = 0
state[mid+1,mid+1] = 1
state[mid+2,mid-1] = 1
return state .> 0
end

function update!(nextstate::BitArray{2},state::BitArray{2})
dimx,dimy = size(state)
for i = 2:dimx-1
for j = 2:dimy -1
count = ifelse(state[i,j],-1,0)
for offi in [-1,0,1]
for offj in [-1,0,1]
if state[i+offi,j+offj]
count += 1
end
end
end

if state[i,j]
if count <2 || count > 3
nextstate[i,j] = false
end
else
if count == 3
nextstate[i,j] = true
end
end
end
end
end

function run(state,nsteps)
nextstate = deepcopy(state)
for i=1:nsteps
update!(nextstate,state)
state  = deepcopy(nextstate)
end
return state
end

nsteps = 1000
n = 100
state = acorn(n)
@time s = run(state,nsteps)
plt.pcolormesh(s)
plt.grid("on")


Output:

  1.746477 seconds (38.55 M allocations: 4.015 GiB, 7.56% gc time)


I see that there is a lot of memory allocation going on. But I am not sure about the best way to avoid this. Ultimately (unless I come up with a completely different algorithm), I do need two matrices (state array and its copy).

What do I have to change in my Julia code to reach similar or better performance than my Python code?

• Did you profile the code to see where are the CPU cycles consumed? Not sure but deepcopy looks like is copy memory in deep and is a candidate for optimize Nov 18, 2020 at 11:29
• The deepcopydoes not seem to be the problem. If I remove it from the inner loop (i.e., I basically repeat the first iteration nsteps times), the timing doesn't change much: 1.649412 seconds (38.42 M allocations: 4.007 GiB, 8.31% gc time). I will have to look into profiling tools, I'm not familiar with them yet. Nov 18, 2020 at 11:40
• Okay it seems that for offi in [-1,0,1] is the bottleneck. I can obviously replace it with for i=-1:1 which speeds up things significantly. 0.202571 seconds (7.01 k allocations: 1.772 MiB) Nov 18, 2020 at 11:45
• n // 2 should become n ÷ 2, I believe. Int(n/2) fails for odd numbers. Dec 19, 2020 at 13:59

Here is my "optimized" version of your code, I'm going down to about 32 7 ms

I use @btime from BenchmarkTools to average multiple trys and avoid counting the compilation time

I made 4 5 changes :

• switch from BitArray to Array{Int} (I don't know why it's faster though)
• update state with .= instead of a copy
• unroll the short loops
• update nextstate with a branchless ifelse
• (EDIT) @inbounds bypasses checking if indices are in the bounds of the array (only use it if you're sure the index isn't greater than the size of the array)

I also changed count to counter since count is a built-in function.

using PyPlot, BenchmarkTools

function acorn(n)
mid = Int(n/2)+1
state = zeros(Int, n, n) # <- specify Int (will be Float64 if not specified)
state[mid,mid-2:mid+4] .= 1
state[mid,mid] = 0
state[mid,mid+1] = 0
state[mid+1,mid+1] = 1
state[mid+2,mid-1] = 1

# I found that Array{Int} is faster than BitArray. Array{Bool} is slowest
return state
end

function update!(nextstate,state)
dimx,dimy = size(state)
for j in 2:dimy-1
for i in 2:dimx-1
# unroll loops
@inbounds counter = state[i-1,j-1] + state[i-1,j] + state[i-1,j+1] +
state[i  ,j-1]                + state[i  ,j+1] +
state[i+1,j-1] + state[i+1,j] + state[i+1,j+1]

# branchless if statement
@inbounds nextstate[i,j] = ifelse(state[i,j] == 1, 1 < counter < 4, counter == 3)

end
end
end

function run(state,nsteps)
nextstate = copy(state)
state = copy(state)

for i in 1:nsteps
update!(nextstate,state)

# each value from nextstate is copied to state, instead of copying the whole array
state .= nextstate
end
return state
end

nsteps = 1000
n = 100
state = acorn(n)
s = @btime run(state,nsteps)
plt.pcolormesh(s)
plt.grid("on")
$$$$

• Thank you! This was the kind of performance I expected but couldn't achieve. Int vs BitArray is strange though. Thanks! Nov 20, 2020 at 8:32
• The downside of BitArrays is that they don't play well with vectorization since computers can only address bytes at a time. This means that re-ordering getindex and setindex` can have really weird results. Nov 21, 2020 at 5:12