1d diffusion in Julia (0.6) and Python (3.5)

test_diffusion.jl

using BenchmarkTools

# jmax iterations of a Laplacian smoother

function smooth!(x,jmax)
    imax = length(x)

    for j = 1:jmax
        xm = x[1]

        for i = 2:imax-1
            xf = 0.5 * x[i] + 0.25 * (x[i+1] + xm)
            xm = x[i]
            x[i] = xf
        end
    end
end

function smooth_v2!(x,jmax)
    imax = length(x)

    xf = similar(x)

    for j = 1:jmax
        for i = 2:imax-1
            xf[i] = 0.5 * x[i] + 0.25 * (x[i+1] + x[i-1])
        end

        # view just gets a reference of the data, it does not copy it
        x[2:imax-1] = view(xf,2:imax-1)
    end

end



function smooth_vec!(x,jmax)
    imax = length(x)

    for j = 1:jmax
        x[2:end-1] = 0.5 * x[2:end-1] + 0.25 * (x[3:end] + x[1:end-2])
    end

end

imax = 1000
x0 = zeros(imax); x0[1] = 1

x = copy(x0)
smooth!(x,100)
@show x[10]

x = copy(x0)
smooth_v2!(x,100)
@show x[10]

x = copy(x0)
smooth_vec!(x,100)
@show x[10]


x = copy(x0)
@btime smooth!(x,100)
#  161.080 μs (0 allocations: 0 bytes)

x = copy(x0)
@btime smooth_v2!(x,100)
#  363.641 μs (101 allocations: 12.63 KiB)

x = copy(x0)
@btime smooth_vec!(x,100)
#  1.174 ms (700 allocations: 5.43 MiB)


# julia: 0.6.2

test_diffusion.py

from __future__ import print_function
import numpy as np
import time


# jmax iterations of a Laplacian smoother

def smooth(x,jmax):
    imax = len(x)

    for j in range(jmax):
        xm = x[0]

        for i in range(1,imax-1):
            xf = 0.5 * x[i] + 0.25 * (x[i+1] + xm)
            xm = x[i]
            x[i] = xf



def smooth_vec(x,jmax):
    imax = len(x)

    for j in range(jmax):
        x[1:imax-1] = 0.5 * x[1:imax-1] + 0.25 * (x[2:imax] + x[0:imax-2])


imax = 1000
x0 = np.zeros((imax,)); x0[0] = 1

x = x0.copy()
start = time.time()
smooth(x,100)
end = time.time()
print("x[9] ", x[9])
print("time (ms) ", (end - start)*1000)
# time (ms)  124.96066093444824


x = x0.copy()
start = time.time()
smooth_vec(x,100)
end = time.time()
print("x[9] ", x[9])
print("time (ms) ", (end - start)*1000)
# time (ms)  1.577615737915039


# python: 3.5.2
# numpy: 1.14.2

Hi, Alexander. I Numba-fied the Pure python version of the above and ran timings using %timeit (friendly wrapper for the timeit module in Python) in a jupyter notebook. I also ran the Julia version above. I'm using Python 2.7.4, Numpy 1.14.5, Numba 0.39, and Julia 0.6.1. Note that the Python packages come from Anaconda and so are built with Intel MKL, while the Julia package is from conda-forge and built with OpenBLAS. Note sure if this makes a difference, but worth noting. Code to generate & test the Numba version is below; note that JIT compilation makes the first execution slow, which timeit ignores:

import numba

smooth_nb = numba.jit(smooth)

x = x0.copy()
%timeit smooth_nb(x, 100)

Version	Time (μs)
Julia	125.809
Pure Python	61400
Numpy	659
Numba	108

Also: you might be interested in the stencil decorator, though I have never used it myself.