6
\$\begingroup\$

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

enter image description here

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)

enter image description here

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?

\$\endgroup\$
4
  • \$\begingroup\$ 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 \$\endgroup\$ – camp0 Nov 18 '20 at 11:29
  • \$\begingroup\$ 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. \$\endgroup\$ – Thomas Nov 18 '20 at 11:40
  • \$\begingroup\$ 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) \$\endgroup\$ – Thomas Nov 18 '20 at 11:45
  • 1
    \$\begingroup\$ n // 2 should become n ÷ 2, I believe. Int(n/2) fails for odd numbers. \$\endgroup\$ – phipsgabler Dec 19 '20 at 13:59
5
\$\begingroup\$

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")
```
\$\endgroup\$
2
  • 1
    \$\begingroup\$ Thank you! This was the kind of performance I expected but couldn't achieve. Int vs BitArray is strange though. Thanks! \$\endgroup\$ – Thomas Nov 20 '20 at 8:32
  • \$\begingroup\$ 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. \$\endgroup\$ – Oscar Smith Nov 21 '20 at 5:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.