# Nbabel nbody integrator speed up

**URL:** <https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712>\
**Category:** Performance\
**Tags:** question, review\
**Created:** [December 12, 2020, 11:26am UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712 "2020-12-12T11:26:19Z")\
**Posts on this page:** 20\
**Page:** 3

<div class="post-metadata">

**Author:** ![paugier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paugier/32/21046_2.png) [@paugier](https://discourse.julialang.org/u/paugier)\
**Post date:** [January 17, 2021, 9:57pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/41 "2021-01-17T21:57:36Z")

</div>

Oh I had a problem to install the nightly rust with rustup. Because it does not seem compatible with rust installed with apt and I’d like to keep this rust for other stuff. And then I was too lazy to do it in a docker 🙂

---

<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:** [January 17, 2021, 10:53pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/42 "2021-01-17T22:53:31Z")

</div>

The greatest difference between the “original” Julia version and the others is that the loop that computed accelerations was not being performed in the same way than in all the other versions. In [this post](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/10) I fixed that. In that version I also use a Julia package called `StaticArrays`, which is useful when one is dealing with many small arrays and want them to be stack-allocated. Not really comparable, but that is just using a package, like one has to use `numpy` properly. The rest is just the natural way to write code in Julia. Of course each one, being more familiar to one or other language will find some things more natural than others. At the final version I don’t use the `StaticArrays` package and manually define the `4-field` struct, but the natural way to go would be to use it.

At the end, the function that computes the accelerations (and take most of the time), is very easy to read, and does not contain anything fancy, I think:

```julia
function compute_acceleration!(positions, masses, accelerations)
    N = length(positions)

    for i in eachindex(accelerations)
        accelerations[i] = zero(eltype(accelerations))
    end

    for i = 1:N - 1
        for j = i + 1:N
            dr = positions[i] - positions[j]
            rinv3 = 1 / norm(dr)^3
            accelerations[i] -= masses[i] * rinv3 * dr
            accelerations[j] += masses[j] * rinv3 * dr
        end
    end
    return nothing
end

```

The “strange” things there are maybe the `eachindex` which can be replaced by `1:length()`, and the `zeros(eltype())`, which could be replaced as well by something syntactically simpler. Those are not for performance, but for generality. That same function can be used for “points” of any dimension, for example, if the appropriate function defining what is the `zero` of such point is defined. Actually it can be used to any kind of `object` (the term does not really apply here) for which the programmer defines what are the zero, the norm, subtraction and the multiplication by a scalar.

---

<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:** [January 18, 2021, 3:57am UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/43 "2021-01-18T03:57:59Z")

</div>

This isn’t an example of advanced high-performance Julia code, except for the `@simd` macros and knowing to use length 4 rather than 3. But it reads pretty clean and simple to me (I mean it’s good code!), but just functions and for-loops. This is also the kind of simple code you should be able to get running fast with numpy as there isn’t much different that the compiler would do.

More advanced high-performance julia you will see much more use of the type system to pass more information to the compiler, maybe `@generated` functions and things like the LoopVectorization.jl `@avx` macro. But as someone who also writes high performance C++, julia is much easier.

---

<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:** [January 18, 2021, 12:59pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/44 "2021-01-18T12:59:30Z")

</div>

Concerning the benchmarks, as for every of these comparisons, one has to be careful in not comparing things improperly. It does not make much sense that a C++ implementation of this is 40% slower than a Fortran implementation of the same thing, also it does not make much sense that any Julia implementation is faster than both.

I have myself fallen into that trap recently, thinking that a Julia code was faster than my own Fortran implementation of the same thing, and an advice of a advanced _Julia_ user taught me that I was simply using the wrong compiler flags for the compilation of the Fortran code. With that fixed the codes ran at the same speed.

There is not much than can be improved in those implementations without going to parallel or cuda versions, thus it does not make much sense that any of the alternatives (Julia, C++, Fortran, Java) are much slower than the other.

It is also natural that the experienced Python developer that implemented the Python version got a better implementation than the sub-optimal implementations of the other languages. I could not implement the same thing in Python as of today.

That being said, of course one would expect that the Python version was slower, and the efficiency that @paugier gets with Transonic shows what a good job they are doing in bringing fast computing to Python. Kudos to them! That certainly puts pressure in all other languages that want to occupy that space! 🙂

---

<div class="post-metadata">

**Author:** ![paugier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paugier/32/21046_2.png) [@paugier](https://discourse.julialang.org/u/paugier)\
**Post date:** [January 18, 2021, 1:30pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/45 "2021-01-18T13:30:42Z")

</div>

> It is also natural that the experienced Python developer that implemented the Python version got a better implementation than the sub-optimal implementations of the other languages. I could not implement the same thing in Python as of today.

@lmiq I didn’t write the C++ and Fortran implementations. These implementations have been written by C++ and Fortran developers and used in the paper mentioned in the README ([GitHub - paugier/nbabel](https://github.com/paugier/nbabel)). They come from the website [https://www.nbabel.org/](https://www.nbabel.org/).

It would be really simple even for me who is not a great C++ dev to improve the C++ version. But this would not be very interesting. Of course, C++ can be very fast.

If a Python implementation using Transonic-Pythran is faster than a C++ implementation, it just means that the C++ implementation is not optimal. Because, behind the scene, a pure C++ implementation (without interaction with CPython) has been written automatically from the Python code by Pythran. It is just much more convenient to write it in Python-Numpy rather than in C++.

Regarding the compiler flags, they are all the same for the different implementations.

---

<div class="post-metadata">

**Author:** ![Skoffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skoffer/32/378_2.png) [@Skoffer](https://discourse.julialang.org/u/Skoffer)\
**Post date:** [January 18, 2021, 2:01pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/46 "2021-01-18T14:01:59Z")

</div>

I vaguely remember that C++ implementation uses N^2 number of calculations for `compute_acceleration` function, while current Julia and python implementations uses only N^2/2 operations. Since `compute_acceleration` dominates in these calculations, it’s not a surprise that C++ and Fortran performs worse compared to Julia and Python versions.

---

<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:** [January 18, 2021, 2:11pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/47 "2021-01-18T14:11:41Z")

</div>

> [@paugier](#):
>
> I didn’t write the C++ and Fortran implementations.

I just took a look at the Fortran implementation, and they use some sort of “large” integer and “float” representations (no idea why). Changing to standard precision appears to have improved a bit the time. Also, the acceleration function there computes the total potential energy as well, while ours (python and Julia versions) do not. So that is “unfair” as well, commenting that line brings the time down another 10%.

Anyway, as I understand the purpose of the repository is to show that Transonic-Pythran can make a Python implementation, which is _expected to be slow_ to run in the same order of the the implementations of languages which are _expected to be fast_. That is clearly demonstrated.

What I do not think is that the remaining differences can be attributed to the languages.

ps. I still cannot run the python version here, I get these errors now:

```julia
conda create -n env_py python pythran transonic pandas
Collecting package metadata (current_repodata.json): done
Solving environment: failed with repodata from current_repodata.json, will retry with next repodata source.
Collecting package metadata (repodata.json): done
Solving environment: failed

PackagesNotFoundError: The following packages are not available from current channels:

  - pythran
  - transonic

Current channels:

  - https://repo.anaconda.com/pkgs/main/linux-64
  - https://repo.anaconda.com/pkgs/main/noarch
  - https://repo.anaconda.com/pkgs/r/linux-64
  - https://repo.anaconda.com/pkgs/r/noarch

To search for alternate channels that may provide the conda package you're
looking for, navigate to

    https://anaconda.org

and use the search bar at the top of the page.

```

---

<div class="post-metadata">

**Author:** ![Skoffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skoffer/32/378_2.png) [@Skoffer](https://discourse.julialang.org/u/Skoffer)\
**Post date:** [January 18, 2021, 2:16pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/48 "2021-01-18T14:16:48Z")

</div>

It could be

```bash
conda create -n nbabel
conda activate nbabel
conda install python=3
conda install -c conda-forge pythran transonic pandas

```

No offense, but I miss Manifest.toml already 🙂

---

<div class="post-metadata">

**Author:** ![Skoffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skoffer/32/378_2.png) [@Skoffer](https://discourse.julialang.org/u/Skoffer)\
**Post date:** [January 18, 2021, 2:24pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/49 "2021-01-18T14:24:56Z")

</div>

Interestingly enough, I’ve got opposite results.

```julia
julia nbabel3_16.jl ../data/input16k 131.16s user 0.54s system 100% cpu 2:11.20 total

python bench.py ../data/input16k 0.2 180.51s user 0.41s system 100% cpu 3:00.36 total

```

so YMMV. I’ve used other julia implementation though: [NBabel with 4th dimension for simd calculations. · GitHub](https://gist.github.com/Arkoniak/d729b119b0bdf218bbf282b847089bca)

---

<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:** [January 18, 2021, 2:40pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/50 "2021-01-18T14:40:42Z")

</div>

Now I get:

```julia
CommandNotFoundError: Your shell has not been properly configured to use 'conda activate'.
To initialize your shell, run...

```

I followed the steps mentioned in the error message, but still nothing… Python is too complicated for me 🤣

---

<div class="post-metadata">

**Author:** ![Skoffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skoffer/32/378_2.png) [@Skoffer](https://discourse.julialang.org/u/Skoffer)\
**Post date:** [January 18, 2021, 2:42pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/51 "2021-01-18T14:42:53Z")

</div>

It should be something like

```julia
source <path to your conda installation>

```

In my case for example

```julia
source ~/miniconda3/etc/profile.d/conda.sh

```

P.S.: I hope admins of this forum wouldn’t be too angry with this python discussion.

---

<div class="post-metadata">

**Author:** ![paugier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paugier/32/21046_2.png) [@paugier](https://discourse.julialang.org/u/paugier)\
**Post date:** [January 18, 2021, 7:36pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/52 "2021-01-18T19:36:08Z")

</div>

I would do once `conda config --add channels conda-forge` and then:

```bash
conda create -n nbabel python pythran transonic pandas
conda activate nbabel

```

We have environment.yaml for conda and pyproject.toml but I didn’t use that for such simple setup.

Don’t worry. I also struggle with Julia installations! In particular, there are many Julia scripts/snippets without instructions for installation in the web.

---

<div class="post-metadata">

**Author:** ![paugier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paugier/32/21046_2.png) [@paugier](https://discourse.julialang.org/u/paugier)\
**Post date:** [January 18, 2021, 7:43pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/53 "2021-01-18T19:43:54Z")

</div>

> [@Skoffer](#):
>
> Interestingly enough, I’ve got opposite results.

I see 4 possibilities:

- I use clang for Pythran and I guess you didn’t use that.

```julia
$ cat ~/.pythranrc                                                                                                                                                          
[pythran]
complex_hook = True
[compiler]
CXX=clang++-10
CC=clang-10
blas=openblas

```

- your Julia version is faster than the one in the repository. If it is the case, PR is really welcome!

- Hardware dependency

- My computer prefers Python and yours prefers Julia just because they know their user.

---

<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:** [January 18, 2021, 8:10pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/54 "2021-01-18T20:10:09Z")

</div>

> [@paugier](#):
>
> My computer prefers Python and yours prefers Julia just because they know their user.

That is why double-blind, placebo trials are necessary.

> [@paugier](#):
>
> `conda activate nbabel`

No luck there, still the same error (the other two installed a bunch of things and did not report anything wrong). This one keeps saying that `CommandNotFoundError: Your shell has not been properly configured to use 'conda activate'`, but using `conda init bash` as indicated does not work (even after restarting the shell, as suggested). ☹

---

<div class="post-metadata">

**Author:** ![paugier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paugier/32/21046_2.png) [@paugier](https://discourse.julialang.org/u/paugier)\
**Post date:** [February 3, 2021, 10:21pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/55 "2021-02-03T22:21:08Z")

</div>

Hello,

I wrote a very simple and naive Julia implementation: [https://github.com/paugier/nbabel/blob/master/julia/nbabel\_naive.jl](https://github.com/paugier/nbabel/blob/master/julia/nbabel_naive.jl)

The goal is to have something as simple and directly comparable to [https://github.com/paugier/nbabel/blob/master/py/bench\_numpy\_highlevel.py](https://github.com/paugier/nbabel/blob/master/py/bench_numpy_highlevel.py)

Can anyone tell me if there is something completely silly/bad in this Julia code?

---

<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 4, 2021, 2:22am UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/56 "2021-02-04T02:22:05Z")

</div>

I don’t see anything clearly bad, although I cannot test it now (the thing we test is if there are type-instabilities, which can make a huge difference, but I don’t see any at first sight).

If the codes are comparable it would depend on what exactly the data structure is when you do this:

```julia
   positions = df.loc[:, ["x", "y", "z"]].values.copy()

```

In the Julia code you are using mutable arrays for representing all data. If the data there is an array of immutable vectors it is likely that performance will be better, and that can be done in Julia with minor changes as well.

(I am not sure if the simd and fast math flags improve anything there, I would leave them out for simplicity, keeping only the inbounds one. The effectiveness of simd operations in particular is very dependent on the data structure, and won’t be great with (3,N) arrays).

---

<div class="post-metadata">

**Author:** ![paugier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paugier/32/21046_2.png) [@paugier](https://discourse.julialang.org/u/paugier)\
**Post date:** [February 4, 2021, 7:54am UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/57 "2021-02-04T07:54:30Z")

</div>

Thanks a lot!

> the thing we test is if there are type-instabilities

Is there a page where I can learn how to do that in Julia?

```julia
positions = df.loc[:, ["x", "y", "z"]].values.copy()

```

creates a C-order contiguous array Nx3 (mutable, equivalent to `Array{Float64}(undef, 3, N)` I guess).

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [February 4, 2021, 8:40am UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/58 "2021-02-04T08:40:02Z")

</div>

The main tool is `@code_warntype`; there’s a discussion of basic usage within the [performance tips](https://docs.julialang.org/en/v1/manual/performance-tips/#man-code-warntype).

---

<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 4, 2021, 9:07am UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/59 "2021-02-04T09:07:55Z")

</div>

> [@paugier](#):
>
> creates a C-order contiguous array Nx3 (mutable, equivalent to `Array{Float64}(undef, 3, N)` I guess).

Is the the same as an `np.array(N,3)`? Why don’t use this?

In the python version you do, for example,

```julia
position0 = positions[index_p0]

```

What is going on there? Is that identical to copying the data to a mutable `np.array(3)` and is it implicit that you are copying over the other dimensions? Why would you do this copy at all in this case instead of referring to the 3 indexes inside the loop?

---

<div class="post-metadata">

**Author:** ![paugier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/paugier/32/21046_2.png) [@paugier](https://discourse.julialang.org/u/paugier)\
**Post date:** [February 4, 2021, 1:24pm UTC](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712/60 "2021-02-04T13:24:33Z")

</div>

I think I should have use (more explicit, but same result):

```python
positions = np.ascontiguousarray(df.loc[:, ["x", "y", "z"]].values)

```

Regarding

```python
position0 = positions[index_p0]

```

it creates a “view” without copy. `positions[index_p0]` is equivalent to `positions[index_p0, :]`. `position0` and `positions` share the same data:

```julia
In [1]: import numpy as np                                                      

In [2]: a = np.arange(6).reshape((3, 2))                                        

In [3]: a                                                                       
Out[3]: 
array([[0, 1],
       [2, 3],
       [4, 5]])

In [4]: b = a[0]                                                                

In [5]: b[:] = 1                                                                

In [6]: a                                                                       
Out[6]: 
array([[1, 1],
       [2, 3],
       [4, 5]])

```

> Why would you do this copy at all in this case instead of referring to the 3 indexes inside the loop?

Shorter, more readable, more general (works for N dim arrays) and no copy. I guess Python-Numpy is designed with the belief that too much loops and indices can decrease readability and maintainability.

But now I realize that the behavior is different in Julia and that one needs the `@view` macro to do that…

```julia
julia> a = reshape(collect(0:5), (2, 3))
2×3 Array{Int64,2}:
 0 2 4
 1 3 5

julia> b = a[:, 1];

julia> b[:] .= 1;

julia> a
2×3 Array{Int64,2}:
 0 2 4
 1 3 5

julia> v = @view a[:, 1];

julia> v[:] .= 1;

julia> a
2×3 Array{Int64,2}:
 1 2 4
 1 3 5

```

Good to know! 🙂

[Previous page](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712.md?page=2)

[Next page](https://discourse.julialang.org/t/nbabel-nbody-integrator-speed-up/51712.md?page=4)