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)
1 Like
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.
Yes, this is clearer thank you. Perfect.
This is what it looks like (@Peter_Adelman 's version), with the structs defined with Parameters, and StaticArrays.
with static arrays
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
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.
1 Like
@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!
1 Like
@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.
1 Like
@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.
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.
1 Like
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.