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.
- If the cell is alive, and exactly two or three neighboring cells are alive, the cell stays alive.
- If the cell is alive and fewer than two or more than three neighboring cells are alive, the cell dies.
- If a cell is dead, and exactly three neighbors are alive, it becomes alive.
- 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?
deepcopy
does not seem to be the problem. If I remove it from the inner loop (i.e., I basically repeat the first iterationnsteps
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. \$\endgroup\$for offi in [-1,0,1]
is the bottleneck. I can obviously replace it withfor i=-1:1
which speeds up things significantly.0.202571 seconds (7.01 k allocations: 1.772 MiB)
\$\endgroup\$n // 2
should becomen ÷ 2
, I believe.Int(n/2)
fails for odd numbers. \$\endgroup\$