# Optimizing the calculation of many distances

**URL:** <https://discourse.julialang.org/t/optimizing-the-calculation-of-many-distances/57343>\
**Category:** Performance\
**Created:** [March 16, 2021, 9:59pm UTC](https://discourse.julialang.org/t/optimizing-the-calculation-of-many-distances/57343 "2021-03-16T21:59:42Z")\
**Posts on this page:** 9\
**Page:** 2

<div class="post-metadata">

**Author:** ![jsbyrnes](https://avatars.discourse-cdn.com/v4/letter/j/e79b87/32.png) [@jsbyrnes](https://discourse.julialang.org/u/jsbyrnes)\
**Post date:** [March 17, 2021, 4:00am UTC](https://discourse.julialang.org/t/optimizing-the-calculation-of-many-distances/57343/22 "2021-03-17T04:00:19Z")

</div>

My goodness, this is faster. Thank you so much. This will be useful, to put it mildly.

(Yes, length(ModSub.x) == length(ModSub.v) at all times)

---

<div class="post-metadata">

**Author:** ![jsbyrnes](https://avatars.discourse-cdn.com/v4/letter/j/e79b87/32.png) [@jsbyrnes](https://discourse.julialang.org/u/jsbyrnes)\
**Post date:** [March 17, 2021, 4:01am UTC](https://discourse.julialang.org/t/optimizing-the-calculation-of-many-distances/57343/23 "2021-03-17T04:01:02Z")

</div>

Maybe I exaggerated some to imply no one is using Julia, but it hasn’t caught on. There are few packages for seismic data, thanks.

---

<div class="post-metadata">

**Author:** ![jsbyrnes](https://avatars.discourse-cdn.com/v4/letter/j/e79b87/32.png) [@jsbyrnes](https://discourse.julialang.org/u/jsbyrnes)\
**Post date:** [March 17, 2021, 4:01am UTC](https://discourse.julialang.org/t/optimizing-the-calculation-of-many-distances/57343/24 "2021-03-17T04:01:54Z")

</div>

Yes, this is clearer thank you. Perfect.

---

<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 17, 2021, 12:02pm UTC](https://discourse.julialang.org/t/optimizing-the-calculation-of-many-distances/57343/25 "2021-03-17T12:02:22Z")

</div>

This is what it looks like (@Peter_Adelman 's version), with the structs defined with Parameters, and StaticArrays.

> **with static arrays**
>
> ```julia
> using Parameters
> using StaticArrays
> 
> @with_kw struct ModSub{N}
> n::Int = N
> x::SVector{N,Float64} = N*rand(SVector{N,Float64}) 
> y::SVector{N,Float64} = 100*rand(SVector{N,Float64}) 
> z::SVector{N,Float64} = 100*rand(SVector{N,Float64}) 
> v::SVector{N,Float64} = rand(SVector{N,Float64}) 
> end
> 
> @with_kw struct ModelConst{N}
> u::ModSub{N} = ModSub{N}() 
> A::ModSub{N} = ModSub{N}()
> t::ModSub{N} = ModSub{N}()
> error::Float64 = 0.01
> BIC::Float64 = 1.0
> fit::Float64 = 1.0
> nfit::Float64 = 1.0
> llh::Float64 = 1.0
> end
> 
> @with_kw struct RayConst
> L::Vector{Float64}
> u::Vector{Float64}
> A::Vector{Float64}
> theta::Vector{Float64}
> phi::Vector{Float64}
> t::Vector{Float64} # dimension 1 to allow mutation
> x::Vector{Float64}
> y::Vector{Float64}
> sta::Float64
> evt::Float64
> pick::Float64
> end
> 
> # Constructor from range r and M
> function RayConst(r,M) 
> lr = length(r)
> RayConst(L = ones(length(r)-1),
> u = ones(length(r)-1),
> A = ones(length(r)-1),
> theta = ones(length(r)-1),
> phi = ones(length(r)-1),
> t = [1.0],
> x = range(0,M,length=lr),
> y = range(0,M,length=lr),
> sta = 1.0,
> evt = 1.0,
> pick = 1.0)
> end
> 
> #
> # Option3 from Peter Adelman
> #
> function v_dist(x::Float64, y::Float64, m::ModSub)
> v_sum = zero(Float64)
> inv_sum = zero(Float64)
> @inbounds for i in 1:length(m.x)
> x_val = (m.x[i] - x)^2
> for j in 1:length(m.y)
> distance = x_val + (m.y[j] - y)^2
> v_sum += m.v[i] / distance
> inv_sum += 1 / distance
> end
> end
> v_sum / inv_sum
> end
> function raytracer!(ray::RayConst, model::ModelConst )
> 
> np = length(ray.x)
> 
> @inbounds for k = 1:np
> u = v_dist(ray.x[k], ray.y[k], model.u)
> 
> A = v_dist(ray.x[k], ray.y[k], model.A)
> 
> t = v_dist(ray.x[k], ray.y[k], model.t)
> if k != np
> ray.u[k] = u
> ray.A[k] = A
> ray.theta[k] = t
> end
> if k != 1
> ray.u[k - 1] += u
> ray.u[k - 1] /= 2
> ray.A[k - 1] += A
> ray.A[k - 1] /= 2
> ray.theta[k - 1] += t
> ray.theta[k - 1] /= 2
> end
> end
> 
> dx = diff(ray.x);
> dy = diff(ray.y);
> 
> @. ray.L = sqrt(dx^2 + dy^2)
> @. ray.phi = atan(dy/dx)
> 
> ray.t[1] = sum(@. ray.L * (ray.u * (1 + ray.A * cos(2*(ray.phi - ray.theta)))))
> 
> return ray
> end
> 
> #
> # Initialize values and run
> #
> 
> n=8
> model = ModelConst{n}( u = ModSub{n}(), 
> A = ModSub{n}(v = 0.1*rand(SVector{n,Float64})),
> t = ModSub{n}(v = π*rand(SVector{n,Float64}) - π) )
> ray = RayConst(1:20,100)
> 
> raytracer!(ray, model)
> 
> using BenchmarkTools
> @btime raytracer!($ray,$model)
> 
> ```

The version not using StaticArrays would be:

> **standard arrays**
>
> ```julia
> using BenchmarkTools
> using Parameters
> using StaticArrays
> 
> @with_kw struct ModSub{N}
> n::Int = N
> x::Vector{Float64} = N*rand(N) 
> y::Vector{Float64} = 100*rand(N) 
> z::Vector{Float64} = 100*rand(N) 
> v::Vector{Float64} = rand(N) 
> end
> 
> @with_kw struct ModelConst{N}
> u::ModSub{N} = ModSub{N}() 
> A::ModSub{N} = ModSub{N}()
> t::ModSub{N} = ModSub{N}()
> error::Float64 = 0.01
> BIC::Float64 = 1.0
> fit::Float64 = 1.0
> nfit::Float64 = 1.0
> llh::Float64 = 1.0
> end
> 
> @with_kw struct RayConst
> L::Vector{Float64}
> u::Vector{Float64}
> A::Vector{Float64}
> theta::Vector{Float64}
> phi::Vector{Float64}
> t::Vector{Float64} # dimension 1 to allow mutation
> x::Vector{Float64}
> y::Vector{Float64}
> sta::Float64
> evt::Float64
> pick::Float64
> end
> 
> # Constructor from range r and M
> function RayConst(r,M) 
> lr = length(r)
> RayConst(L = ones(length(r)-1),
> u = ones(length(r)-1),
> A = ones(length(r)-1),
> theta = ones(length(r)-1),
> phi = ones(length(r)-1),
> t = [1.0],
> x = range(0,M,length=lr),
> y = range(0,M,length=lr),
> sta = 1.0,
> evt = 1.0,
> pick = 1.0)
> end
> 
> #
> # Option3 from Peter Adelman
> #
> function v_dist(x::Float64, y::Float64, m::ModSub)
> v_sum = zero(Float64)
> inv_sum = zero(Float64)
> for i in 1:length(m.x)
> x_val = (m.x[i] - x)^2
> for j in 1:length(m.y)
> distance = x_val + (m.y[j] - y)^2
> v_sum += m.v[i] / distance
> inv_sum += 1 / distance
> end
> end
> v_sum / inv_sum
> end
> 
> function raytracer!(ray::RayConst, model::ModelConst )
> 
> np = length(ray.x)
> 
> @inbounds for k = 1:np
> u = v_dist(ray.x[k], ray.y[k], model.u)
> 
> A = v_dist(ray.x[k], ray.y[k], model.A)
> 
> t = v_dist(ray.x[k], ray.y[k], model.t)
> if k != np
> ray.u[k] = u
> ray.A[k] = A
> ray.theta[k] = t
> end
> if k != 1
> ray.u[k - 1] += u
> ray.u[k - 1] /= 2
> ray.A[k - 1] += A
> ray.A[k - 1] /= 2
> ray.theta[k - 1] += t
> ray.theta[k - 1] /= 2
> end
> end
> 
> dx = diff(ray.x);
> dy = diff(ray.y);
> 
> @. ray.L = sqrt(dx^2 + dy^2)
> @. ray.phi = atan(dy/dx)
> 
> ray.t[1] = sum(@. ray.L * (ray.u * (1 + ray.A * cos(2*(ray.phi - ray.theta)))))
> 
> return ray
> end
> 
> #
> # Initialize values and run
> #
> 
> n=8
> model = ModelConst{n}( u = ModSub{n}(), 
> A = ModSub{n}(v = 0.1*rand(n)),
> t = ModSub{n}(v = π*rand(n) .- π) )
> ray = RayConst(1:20,100)
> 
> raytracer!(ray, model)
> 
> @btime raytracer!($ray,$model)
> 
> ```

The one with StaticArrays seems to be about 20% faster.

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [March 17, 2021, 2:27pm UTC](https://discourse.julialang.org/t/optimizing-the-calculation-of-many-distances/57343/26 "2021-03-17T14:27:17Z")

</div>

@jsbyrnes I used to work for CGG (though not as a seismic engineer). Great to see Julia being used here.  
I can put you in touch with the people at Univ Calgary who did the SeismicJulia packages.  
Looks as if these packages have now been updated as of June 2020, which is a good thing!

---

<div class="post-metadata">

**Author:** ![jsbyrnes](https://avatars.discourse-cdn.com/v4/letter/j/e79b87/32.png) [@jsbyrnes](https://discourse.julialang.org/u/jsbyrnes)\
**Post date:** [March 17, 2021, 3:25pm UTC](https://discourse.julialang.org/t/optimizing-the-calculation-of-many-distances/57343/27 "2021-03-17T15:25:42Z")

</div>

@johnh yea they’ve done a lot and might be good to actually talk to someone; I have a data processing package from the Harvard group when I was doing some ambient noise work, and when I move in the future to some event-based data processing I think the Calgary or Leeds (?) packages might be better. Don’t know yet. But they are at the forefront of this conversion. MATLAB has advantages, but for computationally intensive stuff it’s days in geophysics seem numbered.

But at the moment I’m just moving a monte-carlo routine for tomography from MATLAB to Julia. No data processing is involved. This script is a bit misnamed; it’s really just a line integration routine through a medium defined by scattered points. The actual ray bending will be my next step but the bottleneck will be this function, and already implemented in MATLAB.

@lmiq beautiful, thank you. I appreciate all the effort that’s gone into this. It’s made a big difference.

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [March 17, 2021, 3:39pm UTC](https://discourse.julialang.org/t/optimizing-the-calculation-of-many-distances/57343/28 "2021-03-17T15:39:00Z")

</div>

@jsbyrnes thanks for the reply! I am interested in what you say about tomography. At CGG there was interest in setting up a joint seminar/study day with medical imaging people as there is a lot of commonality in the problems there.

---

<div class="post-metadata">

**Author:** ![Jeff\_Emanuel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeff_emanuel/32/15440_2.png) [@Jeff\_Emanuel](https://discourse.julialang.org/u/Jeff_Emanuel)\
**Post date:** [March 17, 2021, 6:22pm UTC](https://discourse.julialang.org/t/optimizing-the-calculation-of-many-distances/57343/29 "2021-03-17T18:22:15Z")

</div>

Medical imaging has the advantage that they can surround the body with a sources and detectors. Solid earth geophysics also can, but very coarsely. Exploration geophysics is at a distinct disadvantage that one cannot surround the volume of interest. Also raypaths through the body are generally straight. That’s not the case for seismic.

---

<div class="post-metadata">

**Author:** ![jsbyrnes](https://avatars.discourse-cdn.com/v4/letter/j/e79b87/32.png) [@jsbyrnes](https://discourse.julialang.org/u/jsbyrnes)\
**Post date:** [March 18, 2021, 5:11pm UTC](https://discourse.julialang.org/t/optimizing-the-calculation-of-many-distances/57343/30 "2021-03-18T17:11:08Z")

</div>

Yea that’s right it’s not the same; a lot of x-ray stuff can be expressed as a radon transform because the ray path coverage is so good. But its hardly a different problem conceptually.

[Previous page](https://discourse.julialang.org/t/optimizing-the-calculation-of-many-distances/57343.md?page=1)