# Improving an algorithm that compute gps distances

**URL:** <https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213>\
**Category:** Performance\
**Tags:** tullio\
**Created:** [April 26, 2020, 12:46am UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213 "2020-04-26T00:46:03Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [April 26, 2020, 2:31pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/21 "2020-04-26T14:31:10Z")

</div>

> [@Vasily\_Pisarev](#):
>
> On my machine (Ivy Bridge i7-3770, DDR3-1333 RAM), tullio version is 4x slower than @DNF’s with `@avx` .

That’s not great, I wonder why? What versions of Julia and of LoopVectorization was this on? (If you have a minute, a run without loading LoopVectorization would test whether that is the issue, I guess it could be fussy about hardware.)

Thanks for pulling them together. Interesting that Jax uses so little memory, I guess it doesn’t actually materialise `diff_lat` as numpy does.

---

<div class="post-metadata">

**Author:** ![Vasily\_Pisarev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vasily_pisarev/32/7929_2.png) [@Vasily\_Pisarev](https://discourse.julialang.org/u/Vasily_Pisarev)\
**Post date:** [April 26, 2020, 2:44pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/22 "2020-04-26T14:44:21Z")

</div>

Without LoopVectorization, the result is much more decent:

```julia
julia> @btime distances_tullio($a, $b);
551.643 ms (648 allocations: 190.91 MiB)

```

LoopVectorization was v0.7.2

```julia
julia> versioninfo()
Julia Version 1.4.1
Commit 381693d3df* (2020-04-14 17:20 UTC)
Platform Info:
  OS: Linux (x86_64-suse-linux)
  CPU: Intel(R) Core(TM) i7-3770 CPU @ 3.40GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, ivybridge)

```

Jax’s memory footprint is the same as others’, it just works with `Float32` arrays while others with `Float64`.

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [April 26, 2020, 3:01pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/23 "2020-04-26T15:01:07Z")

</div>

Thanks, that’s better. Looking things up, i7-3770 does not have AVX2, perhaps that’s the issue.

> [@Vasily\_Pisarev](#):
>
> Jax’s memory footprint is the same as others’, it just works with `Float32` arrays

I meant that, accounting for 32/64, it’s like the efficient Julia algorithms, and unlike the first one / the numpy version, which make some large intermediate arrays. Perhaps this would not surprise someone who knew more about Jax.

---

<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:** [April 26, 2020, 4:21pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/24 "2020-04-26T16:21:07Z")

</div>

> [@dmolina](#):
>
> No, julia does not know that is a constant, you should explicitly say it:

But it’s written exactly this way in the first post:

> [@David\_Cardozo](#):
>
> const None = [CartesianIndex()]

---

<div class="post-metadata">

**Author:** ![dmolina](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dmolina/32/5246_2.png) [@dmolina](https://discourse.julialang.org/u/dmolina)\
**Post date:** [April 26, 2020, 4:53pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/25 "2020-04-26T16:53:34Z")

</div>

Yes, you are right, my fault.

---

<div class="post-metadata">

**Author:** ![cgarciae](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgarciae/32/14366_2.png) [@cgarciae](https://discourse.julialang.org/u/cgarciae)\
**Post date:** [April 26, 2020, 5:52pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/26 "2020-04-26T17:52:37Z")

</div>

@mcabbott @Vasily_Pisarev I would love to update my test with your new code but I don’t see most of the functions you are using on this thread. Are they in a gist or can you post the final version? Thanks!

---

<div class="post-metadata">

**Author:** ![cgarciae](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgarciae/32/14366_2.png) [@cgarciae](https://discourse.julialang.org/u/cgarciae)\
**Post date:** [April 26, 2020, 6:00pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/27 "2020-04-26T18:00:25Z")

</div>

> [@Tamas\_Papp](#):
>
> It never ceases to amaze me how people expect to use a complex tool without _some_ initial investment. But apparently it is a common phenomenon: the other day I was reading an interview with a trauma surgeon on the rise of injuries from people trying to do some home improvement _with power tools_ during social distancing, with predictable results

This is my initial investment -\_- Also Julia has this on its main page:

> Easy to use  
> Julia has high-level syntax, making it an accessible language for programmers from any background or experience level. Browse the Julia microbenchmarks to get a feel for the language.

So is it a complex tool or is it easy to use?

---

<div class="post-metadata">

**Author:** ![David\_Cardozo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/david_cardozo/32/14344_2.png) [@David\_Cardozo](https://discourse.julialang.org/u/David_Cardozo)\
**Post date:** [April 26, 2020, 6:23pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/28 "2020-04-26T18:23:41Z")

</div>

Thank you all @cgarciae and I were trying to get our hands in Julia for some problem we encounter and making benchmarks against numpy and JAX, hopefully we could learn a lot from this thread.

Thanks

---

<div class="post-metadata">

**Author:** ![cgarciae](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgarciae/32/14366_2.png) [@cgarciae](https://discourse.julialang.org/u/cgarciae)\
**Post date:** [April 26, 2020, 6:25pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/29 "2020-04-26T18:25:33Z")

</div>

> [@dmolina](#):
>
> Jax is great (and [Flax](https://github.com/google-research/flax/tree/prerelease)), but it is not yet to newcomers neither. It is a great tool, but I still prefer Julia.

Hmm. I’d have to argue with this: I took existing Numpy code and just by adding `jax.jit` decorator I got all this speedup, and it also runs on GPU if you install Jax with cuda support. I found it really easy to use.

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [April 26, 2020, 6:30pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/30 "2020-04-26T18:30:49Z")

</div>

> [@cgarciae](#):
>
> but I don’t see most of the functions you are using on this thread. Are they in a gist

Vasily collected lots of them in [this post above](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/20).

BTW, I managed to dig up an old computer with an i5-3427U, on which `distances_threaded_simd` seems to benefit a little from `@avx`, but `distances_tullio` is a disaster. So perhaps this can be narrowed down… LoopVectorization uses CpuId.jl and should ideally detect these things.

---

<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:** [April 26, 2020, 6:35pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/31 "2020-04-26T18:35:31Z")

</div>

> [@cgarciae](#):
>
> just by adding `jax.jit` decorator I got all this speedup

It was claimed upthread that Jax only computed the results lazily, which would mean that the comparison wasn’t very relevant. Can you confirm this?

---

<div class="post-metadata">

**Author:** ![kevbonham](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kevbonham/32/216165_2.png) [@kevbonham](https://discourse.julialang.org/u/kevbonham)\
**Post date:** [April 26, 2020, 6:38pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/32 "2020-04-26T18:38:01Z")

</div>

> [@cgarciae](#):
>
> So is it a complex tool or is it easy to use?

Yes. Both. It’s super easy to use, as evidenced by the fact that I’m a biologist with no CS training and I’m about to be incredibly productive with it. But eeking out every ounce of performance is a different matter. The fact that both simplicity and complexity can exist in the same language is the strength.

---

<div class="post-metadata">

**Author:** ![cgarciae](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgarciae/32/14366_2.png) [@cgarciae](https://discourse.julialang.org/u/cgarciae)\
**Post date:** [April 26, 2020, 6:53pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/34 "2020-04-26T18:53:56Z")

</div>

> [@kevbonham](#):
>
> Yes. Both. It’s super easy to use, as evidenced by the fact that I’m a biologist with no CS training and I’m about to be incredibly productive with it. But eeking out every ounce of performance is a different matter. The fact that both simplicity and complexity can exist in the same language is the strength.

Easy-to-use usually refers to the learning curve. I don’t think Julia is hard, I am just surprised new comers like me get “attacked” for not knowing everything from the begging or reading the whole manual, that just not how you learn these days.

---

<div class="post-metadata">

**Author:** ![Vasily\_Pisarev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vasily_pisarev/32/7929_2.png) [@Vasily\_Pisarev](https://discourse.julialang.org/u/Vasily_Pisarev)\
**Post date:** [April 26, 2020, 7:04pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/35 "2020-04-26T19:04:21Z")

</div>

> [@cgarciae](#):
>
> I would love to update my test with your new code but I don’t see most of the functions you are using on this thread.

The code is under “Function definitions” spoiler in post

> [@Improving an algorithm that compute gps distances](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/20):
>
> On my machine (Ivy Bridge i7-3770, DDR3-1333 RAM), tullio version is 4x slower than @DNF’s with @avx. Maybe that’s hardware issue. Full results: 
> 
> > **Function definitions**
> >
> > using Pkg; pkg"add LoopVectorization https://github.com/mcabbott/Tullio.jl" using Base.Threads, LoopVectorization, Tullio function distances(data1, data2) data1 = deg2rad.(data1) data2 = deg2rad.(data2) lat1 = @view data1[:, 1] lng1 = @view data1[:, 2] lat2 = @view data2[:, 1] lng2 = @view data2[:…

> [@David\_Cardozo](#):
>
> Thank you all @cgarciae and I were trying to get our hands in Julia for some problem we encounter and making benchmarks against numpy and JAX, hopefully we could learn a lot from this thread.

That’s fine. Just next time, please share your benchmarks here when you ask for help and get it. As we see, it benefits everyone that way.

---

<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:** [April 26, 2020, 7:08pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/36 "2020-04-26T19:08:19Z")

</div>

I think the amount of patience depends a bit on how cocksure the newcomer appears to be. I notice that _I_ sometimes get a bit snarky when someone new to the language shows off their–quite understandably–flawed Julia code vs some other language, and proceeds loudly to proclaim how much Julia is behind 😉

And this happens quite frequently.

I think most people who ask for help get very friendly treatment, but “with great confidence comes sharper feedback.”

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [April 26, 2020, 7:09pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/37 "2020-04-26T19:09:28Z")

</div>

I don’t think people are attacking you - but you did put a bunch of posts on social media of benchmarks with some pretty inefficient code…

Julia can require some digging to really tune things, but the end result is typically - things faster then JAX. Sometimes rivaling FORTRAN. So it’s probably best when you run benchmarks but are brand new to using/understanding some of your samples to state that. “I don’t really know julia well but JAX seems fast!”, but that’s not how your tweets read.

There’s a lot of myths the julia community has to continuously battle from the blow-hards in the python community(not everyone is like that - but there’s a lot of reinforcement bias going around). Posts like this just make our lives harder - so that’s probably why you’ve met some friction…

JAX is fast though and yea it is pretty easy to use. But just because code is written one way, JAX, will be doing a lot of optimizations behind your back. Just make sure you compare apples to apples when sharing research publicly or yea people might say “hummmm”…

---

<div class="post-metadata">

**Author:** ![cgarciae](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgarciae/32/14366_2.png) [@cgarciae](https://discourse.julialang.org/u/cgarciae)\
**Post date:** [April 26, 2020, 7:42pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/38 "2020-04-26T19:42:10Z")

</div>

@Vasily_Pisarev You should use `np.asarray` as the end since `np.array` forces a copy of the data. I made these changes and got this numbers using 8 cores:

```julia
distances 1.744 s (40 allocations: 286.18 MiB)
distances_bcast 1.464 s (30 allocations: 95.44 MiB)
distances_threaded 330.340 ms (105 allocations: 190.82 MiB)
distances_threaded_simd 150.763 ms (104 allocations: 190.82 MiB)
dist_np_test 1.413 s (39 allocations: 95.37 MiB)
dist_jax_test 259.303 ms (8 allocations: 320 bytes)

```

code: [final.jl · GitHub](https://gist.github.com/cgarciae/a69fa609f8fcd0aacece92660b5c2315)

Wow! Julia + SIMD is amazing. I am guessing Jax doesn’t use SIMD. As noted on Twitter the allocation numbers for Jax don’t mean anything.

Thanks all! I’ve learned a lot today 😃

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 26, 2020, 9:23pm UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/39 "2020-04-26T21:23:56Z")

</div>

> [@cgarciae](#):
>
> it also runs on GPU if you install Jax with cuda support. I found it really easy to use.

If you make the arrays CuArrays then the broadcasted form will act on the GPU. `]add CuArrays` will even install CUDA for you. I’m curious about the GPU timings.

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [April 27, 2020, 8:31am UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/40 "2020-04-27T08:31:23Z")

</div>

Some further tweaking [over here](https://github.com/mcabbott/Tullio.jl/blob/master/benchmarks/01/distances01.jl) gets `distances_threaded_simd` is down to 65 ms (with Float32, on a 6-core CPU) and `distances_tullio` to 28 ms.

`distances_bcast` does also work on the GPU (31 ms, on an ancient one). As does `distances_tullio`, via KernelAbstractions, but this is very slow right now, not sure why.

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [April 27, 2020, 8:41am UTC](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213/41 "2020-04-27T08:41:43Z")

</div>

In the broadcasted version, won’t ` cos.(lat1) .* cos.(lat2')` cause these cosines to be calculated n^2 instead of n, due to the transpose?

[Previous page](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213.md?page=1)

[Next page](https://discourse.julialang.org/t/improving-an-algorithm-that-compute-gps-distances/38213.md?page=3)