# C++ code much faster than Julia how can I optimize it?

**URL:** <https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868>\
**Category:** Performance\
**Created:** [September 27, 2022, 11:05am UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868 "2022-09-27T11:05:36Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [October 9, 2022, 5:47pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/21 "2022-10-09T17:47:34Z")

</div>

> [@Palli](#):
>
> If it’s just a problem of calling C++

It is not. C++ has different semantics from Julia.  
C++ automatically inserts function calls at key points (e.g., when a function receives an argument by value, or when an object leaves scope).  
Julia does not.  
Thus, by simply defining a few methods in C++, you can automate the memory allocation behavior in any way you like.  
But these function calls must be inserted everywhere manually in Julia to replicate the behavior, meaning the only options in Julia are fully manual memory management, or relying on the GC.

---

<div class="post-metadata">

**Author:** ![Palli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palli/32/3380_2.png) [@Palli](https://discourse.julialang.org/u/Palli)\
**Post date:** [October 9, 2022, 6:11pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/22 "2022-10-09T18:11:28Z")

</div>

For me, 7.7 ms or 0.169365/0.007689 = 22x faster!, so since the original C++ code was “almost 9 times faster” I assume it’s now 2.4x slower, so I advise the thread be solved with your solution! I still think maybe it could be made faster with alternative vector (also in C++), and now I’m thinking what’s wrong with the C++ code… but I’m not going to spend time on that!

I’m not sure maybe this was just an exercise for the new Julia user, he forked some C++ code from a Python project. I think it’s fair game to use code people have helped with on discourse (should give you credit), and if it’s better than Julia clustering packages, hopefully we see it merged into some Julia cluster package.

> [@Elrod](#):
>
> But these function calls must be inserted everywhere manually in Julia to replicate the behavior, meaning the only options in Julia are fully manual memory management, or relying on the GC.

Why I mentioned bittypes. In Julia you can use manual memory management, or if calling C or C++, code using it indirectly, and then no issue with the GC, since it’s bypassed (for vectors of e.g. integers).

I don’t know how Julis’s GC scans stackframes, but I believe it can and does conservatively (at least in some cases, to call C code), so ok if some integers look like pointers.

[https://groups.google.com/g/julia-users/c/JeuJ-QUuLJA/m/xNTgRR31ImcJ](https://groups.google.com/g/julia-users/c/JeuJ-QUuLJA/m/xNTgRR31ImcJ)

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [October 9, 2022, 6:40pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/23 "2022-10-09T18:40:01Z")

</div>

> [@PabloZubieta](#):
>
> `CostCalculator` has abstract fields `Array{Float32}` as opposed `Vector{Float32}`

Ah, that one I didn’t notice. That’s clearly a major source of cycle spending.

So, @Adegasel, whereas in many languages, type-annotation of function arguments is important, that’s not the case in Julia. Unless you have functions with the same name but with different arguments, you need not annotate with types. On the other hand, type annotation in `struct`s is extremely important for performance.

The `Array{Float32}` is an abstract type, it can hold a vector (an `Array{Float32, 1}`), a matrix (an `Array{Float32, 2}`), or a higher dimensional array. All that is known to the compiler when `calc` is compiled is that `cc` holds some kind of arrays of `Float32`. It has to figure out at runtime, every time `cc.cumsum`/`cc.cumsum2` is accessed, the dimension of the array, and figure out the right indexing method for this kind of array.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [October 9, 2022, 6:48pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/24 "2022-10-09T18:48:43Z")

</div>

There’s a bunch of side talk about GCs and other more or less random stuff unrelated to the problem at hand here. As almost always, the solution is just to write completely normal Julia code that follows the performance guide as shown in [C++ code much faster than Julia how can I optimize it? - #20 by pabloferz](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/20). I would suggest the other discussion gets moved to its own thread since it is quite distracting and irrelevant.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 9, 2022, 8:29pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/25 "2022-10-09T20:29:56Z")

</div>

Removing all use of `push!` and nearly all allocations in the inner loops shaves this down to 3ms for me from 5 for the version of @PabloZubieta.

> **KMeans1D**
>
> ```julia
> 
> ```

Edit: somehow this wasn’t the right code. But the one below is.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 9, 2022, 8:41pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/26 "2022-10-09T20:41:06Z")

</div>

And using a `Vector{Int32}` instead of a `Dict` gets it down to 2.1 seconds. So probably \>5x faster than the C++?

> **Kmeans1d no Dict**
>
> ```julia
> module Kmeans1D
> 
> function smawk!(result, rows, cols, allocs, lookup, k)
> length(rows) == 0 && return result
> 
> (; cols_alloc, odd_rows, col_idx_lookup) = allocs[1]
> ## REDUCE
> ncols = 0
> @inbounds for col in cols
> while true
> ncols == 0 && break
> row = rows[ncols]
> if lookup(k, row, col) >= lookup(k, row, cols_alloc[ncols])
> break
> end
> ncols -= 1
> end
> if ncols < length(rows)
> ncols += 1
> cols_alloc[ncols] = col
> end
> end
> _cols = view(cols_alloc, Base.OneTo(ncols))
> 
> # Call recursively on odd-indexed rows
> @inbounds for i in 2:2:length(rows)
> odd_rows[i >> 1] = rows[i]
> end
> 
> smawk!(result, odd_rows, _cols, view(allocs, 2:lastindex(allocs)), lookup, k)
> 
> @inbounds for idx in 1:ncols
> col_idx_lookup[_cols[idx]] = idx
> end
> 
> ## INTERPOLATE
> 
> # Fill-in even-indexed rows
> start = 1
> @inbounds for r in 1:2:length(rows)
> row = rows[r]
> stop = length(_cols) - 1
> if r < (length(rows) - 1)
> stop = col_idx_lookup[result[rows[r + 1]]]
> end
> argmin = _cols[start]
> min = lookup(k, row, argmin)
> for c in start+1:stop+1
> value = lookup(k, row, _cols[c])
> if (c == start) || (value < min)
> argmin = _cols[c]
> min = value
> end
> end
> result[row] = argmin
> start = stop
> end
> 
> return result
> end
> 
> struct CostCalculator
> cumsum::Vector{Float32}
> cumsum2::Vector{Float32}
> 
> function CostCalculator(array, n::Integer)
> cumsum = zeros(n+1)
> cumsum2 = zeros(n+1)
> @inbounds for i in 1:n
> x = array[i]
> cumsum[i+1] = x + cumsum[i]
> cumsum2[i+1] = x * x + cumsum2[i]
> end
> @assert length(cumsum) == length(cumsum2) == n + 1
> return new(cumsum, cumsum2)
> end
> end
> 
> function calc(cc::CostCalculator, i, j)
> if j < i
> return zero(eltype(cc.cumsum))
> end
> 
> @inbounds begin
> mu = (cc.cumsum[j + 1] - cc.cumsum[i]) / (j - i + 1)
> result = cc.cumsum2[j + 1] - cc.cumsum2[i]
> result += (j - i + 1) * (mu * mu)
> result -= (2 * mu) * (cc.cumsum[j + 1] - cc.cumsum[i])
> end
> 
> return result
> end
> 
> struct LookUp
> calculator::CostCalculator
> D::Matrix{Float32}
> end
> 
> function (lu::LookUp)(k, i, j)
> col = min(i, j - 1)
> if col == 0
> col = size(lu.D, 2) + col
> end
> return @inbounds lu.D[k - 1, col] + calc(lu.calculator, j, i)
> end
> 
> function cluster(array, k)
> n = length(array)
> return cluster(array, n, min(k, n))
> end
> 
> function cluster(array, n, k)
> # Sort input array and save info for de-sorting
> sort_idx = sortperm(array)
> undo_sort_lookup = Vector{Int32}(undef, n)
> sorted_array = Vector{Float32}(undef, n)
> 
> @inbounds for i in 1:n
> sorted_array[i] = array[sort_idx[i]]
> undo_sort_lookup[sort_idx[i]] = i
> end
> 
> #Set D and T using dynamic programming algorithm
> cost_calculator = CostCalculator(sorted_array, n)
> D = Matrix{Float32}(undef, k, n)
> T = Matrix{Int32}(undef, k, n)
> lookup = LookUp(cost_calculator, D)
> 
> @inbounds for i in 1:n
> D[1, i] = calc(cost_calculator, 1, i)
> T[1, i] = 1
> end
> 
> row_argmins = Vector{Int32}(undef, n)
> rows = Int32.(1:n)
> cols = Int32.(1:n)
> allocs = NamedTuple{(:odd_rows,:cols_alloc,:col_idx_lookup),Tuple{Vector{Int32},Vector{Int32},Vector{Int32}}}[]
> l = length(rows)
> while true 
> l == 0 && break
> odd_rows = Vector{Int32}(2:2:l)
> cols_alloc = zeros(Int32, l)
> col_idx_lookup = Vector{Int32}(undef, length(rows))
> push!(allocs, (; odd_rows, cols_alloc, col_idx_lookup))
> l ÷= 2
> end
> 
> for k_ in 2:k
> smawk!(row_argmins, rows, cols, allocs, lookup, k_)
> @inbounds for i in 1:n
> argmin = row_argmins[i]
> min = lookup(k_, i, argmin)
> D[k_, i] = min
> T[k_, i] = argmin
> end
> end
> 
> #Extract cluster assignments by backtracking
> centroids = zeros(k)
> sorted_clusters = Vector{Int32}(undef, n)
> t = n + 1
> k_ = k
> n_ = n
> 
> @inbounds while t > 1
> t_ = t
> t = T[k_, n_]
> centroid = 0.0
> for i in t:t_-1
> sorted_clusters[i] = k_
> centroid += (sorted_array[i] - centroid) / (i - t + 1)
> end
> centroids[k_] = centroid
> k_ -= 1
> n_ = t - 1
> end
> 
> clusters = Vector{Int32}(undef, n)
> @inbounds for i in 1:n
> clusters[i] = sorted_clusters[undo_sort_lookup[i]]
> end
> 
> return centroids, clusters
> end
> 
> end # module
> 
> ```

---

<div class="post-metadata">

**Author:** ![jar1](https://avatars.discourse-cdn.com/v4/letter/j/c0e974/32.png) [@jar1](https://discourse.julialang.org/u/jar1)\
**Post date:** [October 9, 2022, 8:56pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/27 "2022-10-09T20:56:46Z")

</div>

I’m seeing lots of `@inbounds for i in 1:n` without any type restriction on the array type which seems unsafe.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 9, 2022, 9:02pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/28 "2022-10-09T21:02:59Z")

</div>

Sure, feel free to add a one-based indexing check at the start, or a conversion step?

But mostly those loops are on `Vector` defined in the code, independent of the original type.

---

<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:** [October 9, 2022, 9:05pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/29 "2022-10-09T21:05:30Z")

</div>

> [@Raf](#):
>
> `stop = col_idx_lookup[result[rows[r + 1]]]`

I don’t see how you can comfortably add `@inbounds` to something like this. I would not recommend `@inbounds` at all, except in tightly controlled situations.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 9, 2022, 9:07pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/30 "2022-10-09T21:07:54Z")

</div>

This mostly isn’t my code, I just stuck with @PabloZubieta mostly putting `@inbounds` everywhere (or how would you see the comparative improvment of removing allocations?). And these acutaly are tightly controlled conditions. All of those objects are allocated and assigned to internally.

Everyone is free to post their own versions here too 😉

(You can remove them all for a time of 2.8 seconds, which is still pretty good. But it seems the algorithm is fine because there are no errors doing this.)

---

<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:** [October 9, 2022, 9:13pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/31 "2022-10-09T21:13:40Z")

</div>

> [@Raf](#):
>
> This mostly isn’t my code

Just commenting on your comment to the comment about the code 😉

With triply nested indexing, you need a pretty steady hand to ensure nothing goes out of bounds. Seems risky _with regular_ `Array`s, let alone general indexable collections.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 9, 2022, 9:14pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/32 "2022-10-09T21:14:55Z")

</div>

They’re not generally indexable collections, they’re basically all internally allocated `Vector` and ranges, except at the outermost level at the start.

---

<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:** [October 9, 2022, 9:40pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/33 "2022-10-09T21:40:16Z")

</div>

I don’t understand why you’re doubling down on this, instead of acknowledging that this is a bad use of `@inbounds` _even with regular Arrays_. It’s a triply nested index! There’s a zillion ways to mess that up. I would never use `@inbounds` in my own code for that.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 9, 2022, 9:48pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/34 "2022-10-09T21:48:35Z")

</div>

Because this is a performance comparison with C++ ?

(and like I said, feel free to post your version, I even timed it and said it was fine, above. I’m just trying to make this run fast and have a little fun here ?)

---

<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:** [October 9, 2022, 10:01pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/35 "2022-10-09T22:01:20Z")

</div>

> [@Raf](#):
>
> feel free to post your version

I don’t have my own version to post, and even if I did, posting it wouldn’t help warn against dangerous practices.

I do think there’s a reckless culture surrounding the use of `@inbounds`, and that’s what I’m reacting to. But since you seem to fundamentally disagree, I’ll give up the back-and-forth.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 9, 2022, 10:07pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/36 "2022-10-09T22:07:19Z")

</div>

I think the original idea with using it is to make this a fair fight with C++, which doesn’t have ~~inbounds~~ bounds checking.

It’s great that you encourage not using inbounds, and I fully agree it should mostly not be used. But in a performance comparison with C++, not using `@inbounds` is giving ourselves an unnecessary handicap. Here, its 25% slower.

---

<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:** [October 9, 2022, 10:10pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/37 "2022-10-09T22:10:49Z")

</div>

> [@Raf](#):
>
> ```julia
> rows = collect(Int32.(1:nrows))
> cols = collect(Int32.(1:ncols))
> 
> ```

There is no need to collect these ranges, as far as I can see, is there?

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 9, 2022, 10:12pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/38 "2022-10-09T22:12:12Z")

</div>

Good catch, that function isn’t actually used anymore, only the in place version where I don’t use `collect`. I’ll edit the code.

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [October 9, 2022, 10:22pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/39 "2022-10-09T22:22:58Z")

</div>

> [@Raf](#):
>
> not using `@inbounds` is giving ourselves an unnecessary handicap.

…unless you use `eachindex`, which you don’t seem to use at all.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [October 9, 2022, 10:27pm UTC](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868/40 "2022-10-09T22:27:28Z")

</div>

This isn’t actally my code, again… I’m just removing the allocations in the code higher up to make it faster, not fixing everything.

In some cases we a looping over a step range, so cant use `eachindex`. But in others, totally `eachindex` can be used, please post a version with that…

[Previous page](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868.md?page=1)

[Next page](https://discourse.julialang.org/t/c-code-much-faster-than-julia-how-can-i-optimize-it/87868.md?page=3)