Thank you for sharing this extremely useful code and project. I am currently using JUDI for full waveform inversion and have encountered an issue that I hope you can help with.

The demo is： using Statistics, Random, LinearAlgebra using JUDI, NLopt, HDF5, SegyIO, PyPlot

Load starting model

if ~isfile("$(JUDI.JUDI_DATA)/overthrust_model_2D.h5")

ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/WaveformInversion.jl/2DFWI/overthrust_model_2D.h5")

end

n, d, o, m0 = read(h5open("model.h5", "r"), "n", "d", "o", "m0") n = reshape(n, 2) n = convert(Vector{Int64}, n) d = reshape(d, 2) d = convert(Vector{Float64}, d) o = reshape(o, 2) o = convert(Vector{Float32}, o) m0 = convert(Matrix{Float32}, m0) model0 = Model((n[1], n[2]), (d[1], d[2]), (o[1], o[2]), m0)

model0 = Model((Int(n[1]), Int(n[2])), (Int(d[1]), Int(d[2])), (Int(o[1]), Int(o[2])), m0)

Replace w/ full path the observed data directory

path_to_data = "/public/home/wanghongzhou/.julia/environments/v1.7/Wiener_Filtered_25Hz/" data_name = "test" # common base name of all shots

Set up out-of-core data container

container = segy_scan(path_to_data, data_name, ["GroupX","GroupY","RecGroupElevation","SourceSurfaceElevation","dt"]) d_obs = judiVector(container)

Set up source

src_geometry = Geometry(container; key="source")

Bound constraints

vmin = 1.4f0 vmax = 6.5f0

Load data and create data vector

if ~isfile("$(JUDI.JUDI_DATA)/overthrust_2D.segy")

ftp_data("ftp://slim.gatech.edu/data/SoftwareRelease/WaveformInversion.jl/2DFWI/overthrust_2D.segy")

end

block = segy_read("marmousi_2D.segy")

d_obs = judiVector(block)

Set up wavelet and source vector

src_geometry = Geometry(container; key = "source", segy_depth_key = "SourceDepth") wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.008f0) # 8 Hz wavelet q = judiVector(src_geometry, wavelet)#

############################################## FWI #################################################

Optimization parameters

maxiter = 1 batchsize = 4 fhistory_SGD = zeros(Float32, maxiter) fTerm = 1e3 gTerm = 1e1

Projection operator for bound constraints

proj(x) = reshape(median([vec(mmin) vec(x) vec(mmax)], dims=2), model0.n)

Define misfit function

function fwi_misfit(model::Model, q::judiVector, d::judiVector; misfit = "L2", compute_gradient = true)

# Set up operators
M = judiModeling(model, q.geometry, d.geometry)
J = judiJacobian(M, q)

# Data residual, function value and gradient
if misfit == "L2"
    r = M*q - d
    f = .5f0*norm(r)^2
    compute_gradient == true && (g = adjoint(J)*r)    # gradient not necessary for line search
elseif misifit == "huber"
    r = M*q - d
    f = eps^2*sqrt(1f0 + dot(r, r)/eps^2) - eps^2
    compute_gradient == true && (g = adjoint(J)*r/sqrt(1f0 + dot(r, r)/eps^2))
else
    throw("Wrong misfit method specified")
end
if compute_gradient == true
    return f, g
else
    return f
end

end

Main loop

for j = 1: maxiter

# select current subset of shots
i = randperm(d_obs.nsrc)[1:batchsize]
f, g = fwi_misfit(model0, q[i], d_obs[i])
println("FWI iteration no: ", j, "; function value: ", f)
fhistory_SGD[j] = f

# linesearch
step = backtracking_linesearch(model0, q[i], d_obs[i], f, g,proj, fwi_misfit; alpha=1f0)

# Update model and bound projection
model0.m = proj(model0.m + reshape(step, model0.n))

if f <= fTerm || norm(g) <= gTerm
    break
end

end

Scanning ... /public/home/admin/.julia/environments/v1.7/Wiener_Filtered_25Hz/test96.segy Scanning ... /public/home/admin/.julia/environments/v1.7/Wiener_Filtered_25Hz/test97.segy Scanning ... /public/home/admin/.julia/environments/v1.7/Wiener_Filtered_25Hz/test98.segy Scanning ... /public/home/admin/.julia/environments/v1.7/Wiener_Filtered_25Hz/test99.segy Building forward operator Operator forward ran in 13.72 s Operator forward ran in 8.98 s Operator forward ran in 8.92 s Operator forward ran in 13.23 s ERROR: LoadError: MethodError: no method matching iterate(::SeisCon) Closest candidates are: iterate(::Union{LinRange, StepRangeLen}) at ~/julia-1.7.2/share/julia/base/range.jl:826 iterate(::Union{LinRange, StepRangeLen}, ::Integer) at ~/julia-1.7.2/share/julia/base/range.jl:826 iterate(::T) where T<:Union{Base.KeySet{<:Any, <:Dict}, Base.ValueIterator{<:Dict}} at ~/julia-1.7.2/share/julia/base/dict.jl:695 ... Stacktrace: [1] copyto!(dest::Vector{Any}, src::SeisCon) @ Base ./abstractarray.jl:890 [2] _collect(cont::UnitRange{Int64}, itr::SeisCon, #unused#::Base.HasEltype, isz::Base.HasLength) @ Base ./array.jl:655 [3] collect(itr::SeisCon) @ Base ./array.jl:649 [4] broadcastable(x::SeisCon) @ Base.Broadcast ./broadcast.jl:704 [5] broadcasted(::Function, ::Matrix{Float32}, ::SeisCon) @ Base.Broadcast ./broadcast.jl:1301 [6] materialize(bc::JUDI.MultiSource) @ JUDI ~/.julia/packages/JUDI/DUfUj/src/TimeModeling/Types/broadcasting.jl:45 [7] -(ms1::judiVector{Float32, Matrix{Float32}}, ms2::judiVector{Float32, SeisCon}) @ JUDI ~/.julia/packages/JUDI/DUfUj/src/TimeModeling/Types/broadcasting.jl:63 [8] fwi_misfit(model::Model, q::judiVector{Float32, Matrix{Float32}}, d::judiVector{Float32, SeisCon}; misfit::String, compute_gradient::Bool) @ Main ~/.julia/environments/v1.7/fwi_2D_obc_sgd.jl:93 [9] fwi_misfit(model::Model, q::judiVector{Float32, Matrix{Float32}}, d::judiVector{Float32, SeisCon}) @ Main ~/.julia/environments/v1.7/fwi_2D_obc_sgd.jl:88 [10] top-level scope @ ~/.julia/environments/v1.7/fwi_2D_obc_sgd.jl:158 in expression starting at /public/home/admin/.julia/environments/v1.7/fwi_2D_obc_sgd.jl:154

It seems that this error is due to the SeisCon type not supporting iteration or broadcasting operations as expected. The issue arises during the calculation of the residual r = M * q - d.

To attempt to solve this issue, I tried converting d to a compatible vector format, but it was not successful. Could you provide some guidance on how to properly handle SeisCon type data so that it can be used with judiVector operations?

Thank you very much for your assistance!(^▽^)

slimgroup / JUDI.jl

MethodError: no method matching iterate(::SeisCon) #247