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

Question

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?

Answer 1

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")
```