[THIS IS WRONG] Optimized 2D Potts model code (derived from https://nbviewer.jupyter.org/gist/genkuroki/d586036fc52dfbc41d5f7e5b25ec0e4a)

potts.jl

using Random: default_rng, seed!, rand!
using RandomNumbers
using StaticArrays
#using Plots
#using BenchmarkTools

β_crit(q) = log(1 + √q)
rand_potts2d(q, m, n=m, rng=default_rng()) = rand(rng, Base.OneTo(Int8(q)), m, n)

function potts2d_ifelse!(q, s, β, niters, rng)
    m, n = size(s)
    prob = @SVector [exp(-β*k) for k in -4:4]
    @fastmath @inbounds @simd for iter in 1:niters
        for j in 1:n 
            for i in 1:m
                let NN = s[ifelse(i == 1, m, i-1), j],
                    SS = s[ifelse(i == m, 1, i+1), j],
                    WW = s[i, ifelse(j == 1, n, j-1)],
                    EE = s[i, ifelse(j == n, 1, j+1)],
                    CT = s[i, j]
                    NV = rand(rng, Base.OneTo(Int8(q)))
                    k  = (  (NN == CT) + (SS == CT) + (WW == CT) + (EE == CT)
                          - (NN == NV) - (SS == NV) - (WW == NV) - (EE == NV) )
                    s[i,j] = ifelse(rand(rng) < prob[k+5], NV, CT)
                end
            end
        end
    end
end

function potts2d_ifelse2!(q, s, β, niters, rng)
    m, n = size(s)
    prob = [exp(-β*k) for k in -4:4]
    k = similar(s)
    r = similar(s)
    p = zeros(m, n)
    @inbounds for iter in 1:niters
        fill!(k, 0)
        rand!(rng, r)
        modp1!(r, q)
        rand!(rng, p)
        # horizontal
        for i in 1:m
            k[i,1] += s[i,1] == s[i,n]
            k[i,1] -= r[i,1] == s[i,n]
            k[i,n] += s[i,n] == s[i,1]
            k[i,n] -= r[i,n] == s[i,1]
        end
        for j in 2:n, i in 1:m
            k[i,j] += s[i,j] == s[i,j-1]
            k[i,j] -= r[i,j] == s[i,j-1]
        end
        for j in 1:n-1, i in 1:m
            k[i,j] += s[i,j] == s[i,j+1]
            k[i,j] -= r[i,j] == s[i,j+1]
        end
        # vertical
        for j in 1:n
            k[1,j] += s[1,j] == s[m,j]
            k[1,j] -= r[1,j] == s[m,j]
            k[m,j] += s[m,j] == s[1,j]
            k[m,j] -= r[m,j] == s[1,j]
        end
        for j in 1:n, i in 2:m
            k[i,j] += s[i,j] == s[i-1,j]
            k[i,j] -= r[i,j] == s[i-1,j]
        end
        for j in 1:n, i in 1:m-1
            k[i,j] += s[i,j] == s[i+1,j]
            k[i,j] -= r[i,j] == s[i+1,j]
        end
        # update states
        for j in 1:n, i in 1:m
            s[i,j] = ifelse(p[i,j] < prob[k[i,j]+5], r[i,j], s[i,j])
        end
    end
end

function modp1!(X, n)
    @inbounds for i in eachindex(X)
        X[i] = mod(X[i], n) + one(eltype(X))
    end
    return X
end


seed!(4649)
q = 3
niters = 10^5
s₀ = rand_potts2d(q, 100)

# compile
potts2d_ifelse!(q, copy(s₀), β_crit(q), 10, default_rng())
potts2d_ifelse2!(q, copy(s₀), β_crit(q), 10, default_rng())

seed!(4649)
@time potts2d_ifelse!(q, copy(s₀), β_crit(q), niters, default_rng())
seed!(4649)
@time potts2d_ifelse2!(q, copy(s₀), β_crit(q), niters, default_rng())

Julia 1.4.2 (Linux, AMD Ryzen 5 2400G)

~/t/pottsmodel $ julia potts.jl
 20.557107 seconds (2 allocations: 9.953 KiB)
 11.231889 seconds (7 allocations: 108.188 KiB)