# General questions from Python user

**URL:** <https://discourse.julialang.org/t/general-questions-from-python-user/55475>\
**Category:** Performance\
**Created:** [February 17, 2021, 5:23pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475 "2021-02-17T17:23:32Z")\
**Posts on this page:** 20\
**Page:** 3

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [February 23, 2021, 1:26pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/41 "2021-02-23T13:26:30Z")

</div>

> [@moeddel](#):
>
> be careful with your timings though. `timeit` does report the mean, whereas `btime` reports the minimum run time.

The way I’ve tried to get a better equivalent for Python is to use `min(timeit.repeat(f(), number=1, repeat=1000))`. So, you get the actual runtime of each trial and then just pick off the minimum.

---

<div class="post-metadata">

**Author:** ![lungben](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lungben/32/12314_2.png) [@lungben](https://discourse.julialang.org/u/lungben)\
**Post date:** [February 23, 2021, 1:26pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/42 "2021-02-23T13:26:42Z")

</div>

You can do the following to take the minimum run time in Python:

```julia
def py_time(f, number=1_000):
    timer = timeit.Timer(f)
    time = min(timer.repeat(number=number))
    return time/number

```

In general, minimum is the better / more stable measure for run time because noise is strictly positive (there are no random effects which make your code run faster).

Edit: @mihalybaci was faster 😉

---

<div class="post-metadata">

**Author:** ![aplavin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aplavin/32/222056_2.png) [@aplavin](https://discourse.julialang.org/u/aplavin)\
**Post date:** [February 23, 2021, 9:00pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/43 "2021-02-23T21:00:00Z")

</div>

It’s easy to give examples when simple Python code (with common numerical libraries, of course) is faster than a direct translation to julia. E.g. see an old topic of mine, [Trig functions very slow](https://discourse.julialang.org/t/trig-functions-very-slow/15335).

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [February 28, 2021, 11:32am UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/44 "2021-02-28T11:32:36Z")

</div>

> [@Eben60](#):
>
> You can write Julia code much like you wrote Python code, and it could be **much** faster (Julia code against native Python loops), or **much** slower (e.g. naive Julia against numpy).

In the course of my weekend procrastinations I contrieved an example to check my statement: some home-made implementation of convolution in both languages.

Make up some signal and kernel data in Python (unimportant part which is not benchmarked):

```julia
import numpy as np

def lorentz(x, gamma, mu=0):
    return (gamma**2/(((x-mu)**2)+gamma**2))/(np.pi*gamma)

def signal(x):
    return lorentz(x, 0.1, -0.3) + lorentz(x, 0.03, 0.1) + lorentz(x, 0.2, 0.25)

def gauss(x, sigma, mu=0.0):
    return np.exp(-(((x-mu)/sigma)**2)/2.0)/(sigma*np.sqrt(2*np.pi))

def kernel(x):
    return gauss(x, 0.01)

vgauss = np.vectorize(gauss)
vlorentz = np.vectorize(lorentz)
vkernel = np.vectorize(kernel)
vsignal = np.vectorize(signal)

k_xmax = 0.05
s_xmax = 0.5

prec = 4
len_s = 1 + 10**prec
q = 2*s_xmax/10**prec
len_k2 = int(k_xmax/q)
len_k = 1 + 2*len_k2

s_x = np.linspace(-s_xmax, s_xmax, len_s)
k_x = np.linspace(-k_xmax, k_xmax, len_k)

s_y = vsignal(s_x)
k_y = vkernel(k_x)

```

Now the same in Julia:

```julia
lorentz(x, γ, μ=0) = (γ^2/((x-μ)^2+γ^2))/(π*γ)
signal(x) = lorentz(x, 0.1, -0.3) + lorentz(x, 0.03, 0.1) + lorentz(x, 0.2, 0.25)

gauss(x, σ, μ=0.0) = exp(-(((x-μ)/σ)^2)/2)/(σ*sqrt(2π))
kernel(x) = gauss(x, 0.01)

s_xmax = 0.5
k_xmax = 0.05

prec = 4
len_s = 1 + 10^prec
q = 2*s_xmax/10^prec
len_k2 = Int(k_xmax/q)
len_k = 1 + 2len_k2

s_x = LinRange(-s_xmax, s_xmax, len_s)
k_x = LinRange(-k_xmax, k_xmax, len_k)

s_y = signal.(s_x)
k_y = kernel.(k_x)

```

Now the convolution in Python:

```julia
def convolve1(s_y, k_y, len_s, len_k, q):
    r = np.zeros_like(s_y)
    k_rev = np.flip(k_y)
    len_k2 = int((len_k - 1) / 2)
    fst = 1 + len_k2
    lst = len_s - len_k2 - 1
    for i in range(fst, lst):
        r[i-1] = q * np.sum(s_y[i-len_k2:i-len_k2+len_k] * k_rev)
    return r

```

and in Julia:

```julia
function convolve1(s_y, k_y, len_s, len_k, q)
    r = zeros(len_s)
    k_rev = reverse(k_y)
    len_k2 = Int((len_k - 1) / 2)
    fst = 1 + len_k2
    lst = len_s - len_k2 - 1
    for i in fst:lst
        r[i] = q * sum(s_y[i-len_k2+1:i-len_k2+len_k] .* k_rev)
    end
    return r
end

```

In fact I of course wrote Julia first, but let’s pretend the code was ported from Python to Julia. As to the code: I can’t guarantee it’s mathematically or otherwise correct, but we don’t need no convolution - just some identical code in both languages, reasonably vectorized in Numpy on the Python side - which it is, I believe.

Let’s test:

```julia
def cnv1():
    return convolve1(s_y, k_y, len_s, len_k, q)

import timeit

tm = np.amin(timeit.repeat(cnv1, repeat=50, number=1))

print("tm = ", tm)
# tm = 0.062673

```

```julia
julia> using BenchmarkTools
julia> @btime convolve1(s_y, k_y, len_s, len_k, q);
# 24.115 ms (18003 allocations: 140.71 MiB)

```

Julia is 2.5 times faster - not bad! But still we see a lot of allocations meaning a room for improvement. As mentioned in one of the posts above, let’s use views instead of slices :

```julia
function convolve2(s_y, k_y, len_s, len_k, q)
    r = zeros(len_s)
    k_rev = reverse(k_y)
    len_k2 = Int((len_k - 1) / 2)
    fst = 1 + len_k2
    lst = len_s - len_k2 - 1
    for i in fst:lst
        r[i] = q * sum(view(s_y, i-len_k2+1:i-len_k2+len_k) .* k_rev)
    end
    return r
end

# julia> @btime convolve2(s_y, k_y, len_s, len_k, q);
# 13.560 ms (9003 allocations: 70.40 MiB)

```

Substantial improvement, still a lot of allocations. Next try:

```julia
function convolve3(s_y, k_y, len_s, len_k, q)
    r = zeros(len_s)
    t = similar(k_y) # pre-allocate a buffer
    k_rev = reverse(k_y)
    len_k2 = Int((len_k - 1) / 2)
    fst = 1 + len_k2
    lst = len_s - len_k2 - 1
    for i in fst:lst
        t .= view(s_y, i-len_k2+1:i-len_k2+len_k) .* k_rev
        r[i] = q * sum(t)
    end
    return r
end

# julia> @btime convolve3(s_y, k_y, len_s, len_k, q);
# 1.860 ms (4 allocations: 94.27 KiB)

```

That looks good! Actually I would wish the compiler was smart enough to do the preallocation on it’s own here.

Probably it can be made still faster by some kind of multithreading, as the `for` loop here would perfectly lend itself to parallel execution. Lets’s try it:  
_ **Edit** : don’t try it - see posts below concerning race conditions_

```julia
function convolve4(s_y, k_y, len_s, len_k, q)
    r = zeros(len_s)
    t = similar(k_y)
    k_rev = reverse(k_y)
    len_k2 = Int((len_k - 1) / 2)
    fst = 1 + len_k2
    lst = len_s - len_k2 - 1
    Threads.@threads for i in fst:lst
        t .= view(s_y, i-len_k2+1:i-len_k2+len_k) .* k_rev
        r[i] = q * sum(t)
    end
    return r
end

# julia> @btime convolve4(s_y, k_y, len_s, len_k, q);
# 4.737 ms (25 allocations: 97.73 KiB)

```

It got worse ☹ Maybe the dataset to small for effective use of parallel processing? Let’s make it 10 times bigger:

`prec = 5`

```julia
julia> @btime convolve3(s_y, k_y, len_s, len_k, q);
  305.176 ms (6 allocations: 937.92 KiB)
julia> @btime convolve4(s_y, k_y, len_s, len_k, q);
  216.561 ms (27 allocations: 941.39 KiB)

```

Here we indeed got some improvement by multithreading. Can it be improved further?

Now we got it much better, let’s try to get it much worse 😉 Let’s feed our first implementation not with arrays of floats, but with arrays of `Any`. It is not as contrieved as some would probably claim: The data could come e.g. from a file or some other source without the type guarantee.

```julia
s_y_any = Array{Any,1}(undef, len_s)
s_y_any .= s_y
k_y_any = Array{Any,1}(undef, len_k)
k_y_any .= k_y

# julia> @btime convolve1(s_y_any, k_y_any, len_s, len_k, q);
# 476.886 ms (13499992 allocations: 347.25 MiB)

```

That’t 8 times slower than Python.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 28, 2021, 1:06pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/45 "2021-02-28T13:06:45Z")

</div>

> [@Eben60](#):
>
> ```julia
> t .= view(s_y, i-len_k2+1:i-len_k2+len_k) .* k_rev
> r[i] = q * sum(t)
> 
> ```

I think you can avoid the temporary vector by using a dot product instead (untested):

```julia
r[i] = (view(s_y, i-len_k2+1:i-len_k2+len_k)' * k_rev) * q

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 28, 2021, 1:12pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/46 "2021-02-28T13:12:11Z")

</div>

Also I don’t think you can use multithreading as long as you use a temporary `t`, since there will be race conditions for writing into `t`.

But try threading with the dot product instead.

_Update:_ I indeed get a nice speedup with threading, if I remove the temporary `t`:

```julia
function convolve5t(s_y, k_y, len_s, len_k, q)
    r = zeros(len_s)
    k_rev = reverse(k_y)
    len_k2 = div((len_k - 1), 2) # use div, not Int(x/2), here
    fst = 1 + len_k2
    lst = len_s - len_k2 - 1
    Threads.@threads for i in fst:lst
        r[i] = q * (view(s_y, i-len_k2+1:i-len_k2+len_k)' * k_rev)
    end
    return r
end

```

```julia
jl> @btime convolve5t($s_y, $k_y, $len_s, $len_k, $q);
  419.700 μs (84 allocations: 93.52 KiB)

```

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [February 28, 2021, 1:59pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/47 "2021-02-28T13:59:14Z")

</div>

> [@Eben60](#):
>
> Actually I would wish the compiler was smart enough to do the preallocation on it’s own here.

Actually there you don’t need that allocation at all:

```julia
function convolve_new(s_y, k_y, len_s, len_k, q)
    r = zeros(len_s)
    k_rev = reverse(k_y)
    len_k2 = Int((len_k - 1) / 2)
    fst = 1 + len_k2
    lst = len_s - len_k2 - 1
    for i in fst:lst
        k = firstindex(k_rev)
        t = 0.
        @inbounds @simd for j in i-len_k2+1:i-len_k2+len_k
          t += s_y[j]*k_rev[k]
          k += 1
        end
        r[i] = q*t
    end
    return r
end

```

And this turns out to be 2 times faster than `coevolve3`:

```julia
julia> @btime convolve3($s_y, $k_y, $len_s, $len_k, $q);
  1.665 ms (4 allocations: 94.27 KiB)

julia> @btime convolve_new($s_y, $k_y, $len_s, $len_k, $q);
  880.245 μs (3 allocations: 86.27 KiB)

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 28, 2021, 2:01pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/48 "2021-02-28T14:01:54Z")

</div>

> [@DNF](#):
>
> I think you can avoid the temporary vector by using a dot product

Yes, just doing

```julia
@views function convolve_dot(s_y, k_y, len_s, len_k, q)
    r = zeros(len_s)
    k_rev = reverse(k_y)
    len_k2 = (len_k - 1) ÷ 2
    fst = 1 + len_k2
    lst = len_s - len_k2 - 1
    for i in fst:lst
        r[i] = q * (s_y[i-len_k2+1:i-len_k2+len_k]' * k_rev)
    end
    return r
end

@btime convolve_dot($s_y, $k_y, $len_s, $len_k, $q);

```

it got about 30x faster than the original `convolve1` function on my laptop.

“Avoid allocations in inner loops” is a pretty simple rule, but it is unfamiliar to most people coming from Python or Matlab or R where allocations are a tertiary concern and vectorization is everything.

(It also wouldn’t be terrible to simply write out the inner loop, though to get peak performance one would need `@simd` or `@avx`. This is especially useful if you are doing something more complicated than a convolution. You don’t always need or want to dig through the standard library for a function that does exactly what you need!)

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [February 28, 2021, 2:05pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/49 "2021-02-28T14:05:02Z")

</div>

> [@stevengj](#):
>
> It also wouldn’t be terrible to simply write out the inner loop, though to get peak performance one would need `@simd`

Maybe that is what I did there. Here it is slightly faster than using the dot product:

```julia
julia> @btime convolve_dot($s_y, $k_y, $len_s, $len_k, $q);
  938.193 μs (3 allocations: 86.27 KiB)

julia> @btime convolve_new($s_y, $k_y, $len_s, $len_k, $q);
  880.281 μs (3 allocations: 86.27 KiB)

```

Edit: And that version, threaded (4 threads) works pretty well:

```julia
function convolve_new_threaded(s_y, k_y, len_s, len_k, q)
    r = zeros(len_s)
    k_rev = reverse(k_y)
    len_k2 = Int((len_k - 1) / 2)
    fst = 1 + len_k2
    lst = len_s - len_k2 - 1
    Threads.@threads for i in fst:lst
        k = firstindex(k_rev)
        t = 0.
        @inbounds @simd for j in i-len_k2+1:i-len_k2+len_k
          t += s_y[j]*k_rev[k]
          k += 1
        end
        r[i] = q*t
    end
    return r
end

```

```julia
julia> Threads.nthreads()
4

julia> @btime convolve_new_threaded($s_y, $k_y, $len_s, $len_k, $q);
  249.114 μs (24 allocations: 89.30 KiB)

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 28, 2021, 2:16pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/50 "2021-02-28T14:16:11Z")

</div>

Also keep in mind that removing that temporary array is not just about performance, but also correctness. The different threads will trample each other accessing and modifying the temporary out of sequence. It can even be modified by one thread while another is trying sum it.

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [February 28, 2021, 6:01pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/51 "2021-02-28T18:01:43Z")

</div>

Yes, that is very true! My `convolve4` actually produces a wrong result.  
Somehow I kept it in mind with the race condition, but then forgot.

The following works correctly:

```julia
function convolve6(s_y, k_y, len_s, len_k, q)
    r = zeros(len_s)
    t0 = similar(k_y)
    t = [t0]
    for i in 2:Threads.nthreads() # nthreads is 4 on my computer
        push!(t, copy(t0))
    end

    k_rev = reverse(k_y)
    len_k2 = Int((len_k - 1) / 2)
    fst = 1 + len_k2
    lst = len_s - len_k2 - 1
    Threads.@threads for i in fst:lst
        ti = t[Threads.threadid()]
        ti .= view(s_y, i-len_k2+1:i-len_k2+len_k) .* k_rev
        r[i] = q * sum(ti)
    end
    return r
end
# julia> @btime convolve6($s_y, $k_y, $len_s, $len_k, $q);
# 503.801 μs (30 allocations: 121.88 KiB)

```

(as compared with 1.86 ms for the non-threaded version` convolve3` on my computer).

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [February 28, 2021, 6:18pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/52 "2021-02-28T18:18:35Z")

</div>

now all results on the same (my) computer:

1. Initial variant: `convolve1` # 24.115 ms (18003 allocations: 140.71 MiB)

2. Views instead of slices: `convolve2` # 13.560 ms (9003 allocations: 70.40 MiB)

3. Preallocated array: `convolve3` # 1.860 ms (4 allocations: 94.27 KiB)

4. Preallocated array with `@threads`: `convolve6` # 503.801 μs (30 allocations: 121.88 KiB)

5. Using dot-product with `@threads`: convolve5t # 277.300 μs (24 allocations: 89.73 KiB)

6. Write out the inner loop using `@simd`: `convolve_new` # 982.099 μs (3 allocations: 86.27 KiB)

7. The same plus threading: `convolve_new_threaded` # 256.500 μs (24 allocations: 89.73 KiB)

@DNF and @lmiq , thank you a lot for your suggestions!

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [February 28, 2021, 6:22pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/53 "2021-02-28T18:22:55Z")

</div>

> [@Eben60](#):
>
> thank you a lot for your suggestions!

I am curious how a `numba` version of the python code would perform there.

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [February 28, 2021, 6:26pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/54 "2021-02-28T18:26:17Z")

</div>

or `numpy.convolve(a , v , mode)`

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [February 28, 2021, 7:03pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/55 "2021-02-28T19:03:30Z")

</div>

> [@Eben60](#):
>
> numpy.convolve(a , v , mode)

Well, that is a curiosity (it is probably fast) but doesn’t have the same meaning of comparison relative to a code one can write directly.

My interest is more in these python JIT accelerators, and how they behave relative to Julia code. I will try that, but I may make naive mistakes in python and thus produce meaningless comparisons.

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [February 28, 2021, 7:46pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/56 "2021-02-28T19:46:59Z")

</div>

Well, actually I didn’t expect that python would be that much slower that Julia in this case, as most of the computation is actually vectorized in numpy. But it is possible that it is 10000 loops of pure python is what takes the biggest part of these 60 ms - in this case numba acceleration may help.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [February 28, 2021, 10:20pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/57 "2021-02-28T22:20:05Z")

</div>

Adding `import numba` and `@numba.njit` before every function definition takes the time  
of the python version from `51ms` to `16ms` here (I take no responsibility for that result, just added the flags - I have no experience with that).

> **Code**
>
> ```julia
> import numpy as np
> import numba
> 
> @numba.njit
> def lorentz(x, gamma, mu=0):
> return (gamma**2/(((x-mu)**2)+gamma**2))/(np.pi*gamma)
> 
> @numba.njit
> def signal(x):
> return lorentz(x, 0.1, -0.3) + lorentz(x, 0.03, 0.1) + lorentz(x, 0.2, 0.25)
> 
> @numba.njit
> def gauss(x, sigma, mu=0.0):
> return np.exp(-(((x-mu)/sigma)**2)/2.0)/(sigma*np.sqrt(2*np.pi))
> 
> @numba.njit
> def kernel(x):
> return gauss(x, 0.01)
> 
> vgauss = np.vectorize(gauss)
> vlorentz = np.vectorize(lorentz)
> vkernel = np.vectorize(kernel)
> vsignal = np.vectorize(signal)
> 
> k_xmax = 0.05
> s_xmax = 0.5
> 
> prec = 4
> len_s = 1 + 10**prec
> q = 2*s_xmax/10**prec
> len_k2 = int(k_xmax/q)
> len_k = 1 + 2*len_k2
> 
> s_x = np.linspace(-s_xmax, s_xmax, len_s)
> k_x = np.linspace(-k_xmax, k_xmax, len_k)
> 
> s_y = vsignal(s_x)
> k_y = vkernel(k_x)
> 
> @numba.njit
> def convolve1(s_y, k_y, len_s, len_k, q):
> r = np.zeros_like(s_y)
> k_rev = np.flip(k_y)
> len_k2 = int((len_k - 1) / 2)
> fst = 1 + len_k2
> lst = len_s - len_k2 - 1
> for i in range(fst, lst):
> r[i-1] = q * np.sum(s_y[i-len_k2:i-len_k2+len_k] * k_rev)
> return r
> 
> @numba.njit
> def cnv1():
> return convolve1(s_y, k_y, len_s, len_k, q)
> 
> import timeit
> 
> tm = np.amin(timeit.repeat(cnv1, repeat=50, number=1))
> 
> print("tm = ", tm)
> 
> ```

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [March 6, 2021, 8:41pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/58 "2021-03-06T20:41:39Z")

</div>

If I take your code with inner loop `convolve_new` and remove `@inbounds @simd` , then I get

```julia
# julia> @btime convolve_11($s_y, $k_y, $len_s, $len_k, $q);
# 10.022 ms (3 allocations: 86.27 KiB)

```

which is similar to your numba-JIT python result, so my guess is that numpy-broadcasted multiplication of arrays probably doesn’t use SIMD-vectorization.

Actually I’ve written a similar code myself in the beginning, but with a small but important difference (see below), got a poor result and didn’t put it into discussion especially as I wanted to compare possibly similar Python and Julia codes.

```julia
function convolve_12(s_y, k_y, len_s, len_k, q)
 # ...........#
    for i in fst:lst
        k = firstindex(k_rev)
        # t = 0. # not using it, but:
        for j in i-len_k2+1:i-len_k2+len_k
          r[i] += s_y[j]*k_rev[k] # summation directly into r elements
          k += 1
        end
        r[i] *= q
    end
    return r
end

# julia> @btime convolve_12($s_y, $k_y, $len_s, $len_k, $q);
# 22.445 ms (3 allocations: 86.27 KiB)

```

By putting ´@inbounds @simd´ back (actually ´@inbounds´ saving is minimal here) we win back about the same 10 ms as in your code, just the difference is not from 10 ms to 1 ms, but from 20 ms to 10ms

```julia
        @inbounds @simd for j in i-len_k2+1:i-len_k2+len_k

# julia> @btime convolve_13($s_y, $k_y, $len_s, $len_k, $q);
# 9.928 ms (3 allocations: 86.27 KiB)

```

Two questions:

- Why is it that bad to perform summation directly into array elements?
- Is there anything on this point in the Performance Tips

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [March 6, 2021, 9:55pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/59 "2021-03-06T21:55:17Z")

</div>

> [@Eben60](#):
>
> Why is it that bad to perform summation directly into array elements?

Because then it has to mutate an element in the heap at every iteration of that loop. Instead, with the temporary variable, the value is accumulated in an immutable stack (or cache) variable, which is much faster. The compiler could eventually figure that out and do that automatically, but that is less predictable, it depends on the context of the code. I am guessing that it is not doing that because the index of `r` is part of the counter of the loop where it is modified. But I am just guessing.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 8, 2021, 10:49pm UTC](https://discourse.julialang.org/t/general-questions-from-python-user/55475/60 "2021-03-08T22:49:39Z")

</div>

> [@lmiq](#):
>
> Instead, with the temporary variable, the value is accumulated in an immutable stack (or cache) variable, which is much faster.

(Probably in a CPU register, which is even faster.)

[Previous page](https://discourse.julialang.org/t/general-questions-from-python-user/55475.md?page=2)