FWI with user provided misfit function

[asidosaputra]

Hi, I have tried FWI with user-provided misfit function
I follow the tutorial in this link https://slimgroup.github.io/JUDI.jl/dev/tutorials/05_custom_misfit/

using SegyIO, HDF5, PyPlot, JUDI, Random, LinearAlgebra, Printf, SlimOptim, Statistics

sx(d::judiVector) = d.geometry.xloc[1][1]

m0, n, d, o = read(h5open("/media/swina/WORK_DISK_1/Julia/FWI-Julia/scripts/overthrust_model.h5","r"),"m0","n","d","o"); #title("Starting model")
model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)


block = segy_read("/media/swina/WORK_DISK_1/Julia/FWI-Julia/scripts/overthrust_shot_records.segy");
d_obs = judiVector(block);

src_geometry = Geometry(block; key="source");
src_data = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.008f0);
q = judiVector(src_geometry, src_data);

function mse(dsyn, dobs)
    fval = .5f0 * norm(dsyn - dobs)^2
    adjoint_source = syn - dobs
    return fval, adjoint_source
end

function studentst(dsyn, dobs; k=2)
    fval = sum(.5 * (k + 1) * log.(1 .+ (dsyn - dobs).^2 ./ k))
    adjoint_source = (k + 1) .* (dsyn - dobs) ./ (k .+ (dsyn - dobs).^2)
    return fval, adjoint_source
end

d_outlier = deepcopy(d_obs)

## Add outliers to the data
for s=1:d_outlier.nsrc
    # wrongly scale 10 traces
    nrec = d_outlier[s].geometry.nrec[1]
    inds = rand(1:nrec, 10)
    d_outlier.data[s][:, inds] .*= 20
end

g_const = 0
function objective_function(x, d_obs, misfit=mse)
    model0.m .= reshape(x,model0.n);

    # fwi function value and gradient
    i = randperm(d_obs.nsrc)[1:batchsize]
    fval, grad = fwi_objective(model0, q[i], d_obs[i]; misfit=misfit)
    # Normalize for nicer convergence
    g = grad.data
    # Mute water
    g[:, 1:20] .= 0f0
    g_const == 0 && (global g_const = 1/norm(g, Inf))
    return fval, g_const .* g
end

# Bound constraints based on initial velocity
mmin = .95f0 * minimum(m0) * ones(Float32, n...)
mmax = 1.1f0 * maximum(m0) * ones(Float32, n...)

# Fix water layer
mmin[:, 1:20] .= m0[1,1]
mmax[:, 1:20] .= m0[1,1]

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

batchsize = 8
niter = 20

# Setup SPG
options = spg_options(verbose=3, maxIter=niter, memory=3)

# Compare l2 with students t on ideal data
ϕmse_ideal = x->objective_function(x, d_obs)
ϕst_ideal = x->objective_function(x, d_obs, studentst)

# Compare l2 with students t on the data with outliers
ϕmse_out = x->objective_function(x, d_outlier)
ϕst_out = x->objective_function(x, d_outlier, studentst)

# Perform the inversion
g_const = 0
solmse_ideal = spg(ϕmse_ideal, m0, proj, options)
g_const = 0
solst_ideal = spg(ϕst_ideal, m0, proj, options)
g_const = 0
solmse_out = spg(ϕmse_out, m0, proj, options)
g_const = 0
solst_out = spg(ϕst_out, m0, proj, options)

After running all code, i have error below:

julia> include("scripts/02_res_def.jl")
┌ Warning: Fixed length trace flag set in stream: IOBuffer(data=UInt8[...], readable=true, writable=false, seekable=true, append=false, size=7076688, maxsize=Inf, ptr=3601, mark=-1)
└ @ SegyIO ~/.julia/packages/SegyIO/qkvUT/src/read/read_file.jl:26
ERROR: LoadError: MethodError: no method matching multi_src_fg(::Model, ::judiVector{Float32, Matrix{Float32}}, ::judiVector{Float32, Matrix{Float32}}, ::Nothing, ::JUDIOptions, ::Bool, ::Bool, ::typeof(mse))
Closest candidates are:
  multi_src_fg(::Model, ::judiVector, ::judiVector, ::Any, ::JUDIOptions, ::Bool, ::Bool) at ~/.julia/packages/JUDI/RHhnq/src/TimeModeling/Modeling/misfit_fg.jl:4
Stacktrace:
  [1] run_and_reduce(func::Function, #unused#::Nothing, nsrc::Int64, arg_func::JUDI.var"#209#210"{JUDIOptions, Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:nlind, :lin, :misfit), Tuple{Bool, Bool, typeof(mse)}}}, Model, judiVector{Float32, Matrix{Float32}}, judiVector{Float32, Matrix{Float32}}, Nothing})
    @ JUDI ~/.julia/packages/JUDI/RHhnq/src/TimeModeling/Modeling/propagation.jl:37
  [2] multi_src_fg!(G::PhysicalParameter{Float32}, model::Model, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}, dm::Nothing; options::JUDIOptions, kw::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:nlind, :lin, :misfit), Tuple{Bool, Bool, typeof(mse)}}})
    @ JUDI ~/.julia/packages/JUDI/RHhnq/src/TimeModeling/Modeling/propagation.jl:86
  [3] #multi_exp_fg!#192
    @ ~/.julia/packages/JUDI/RHhnq/src/TimeModeling/Modeling/misfit_fg.jl:155 [inlined]
  [4] fwi_objective!(G::PhysicalParameter{Float32}, model::Model, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}; options::JUDIOptions, kw::Base.Pairs{Symbol, typeof(mse), Tuple{Symbol}, NamedTuple{(:misfit,), Tuple{typeof(mse)}}})
    @ JUDI ~/.julia/packages/JUDI/RHhnq/src/TimeModeling/Modeling/misfit_fg.jl:136
  [5] fwi_objective(model::Model, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}; options::JUDIOptions, kw::Base.Pairs{Symbol, typeof(mse), Tuple{Symbol}, NamedTuple{(:misfit,), Tuple{typeof(mse)}}})
    @ JUDI ~/.julia/packages/JUDI/RHhnq/src/TimeModeling/Modeling/misfit_fg.jl:101
  [6] objective_function(x::Matrix{Float32}, d_obs::judiVector{Float32, Matrix{Float32}}, misfit::Function)
    @ Main /media/swina/WORK_DISK_1/Julia/FWI-Julia/scripts/02_res_def.jl:44
  [7] objective_function(x::Matrix{Float32}, d_obs::judiVector{Float32, Matrix{Float32}})
    @ Main /media/swina/WORK_DISK_1/Julia/FWI-Julia/scripts/02_res_def.jl:40
  [8] (::var"#3#4")(x::Matrix{Float32})
    @ Main /media/swina/WORK_DISK_1/Julia/FWI-Julia/scripts/02_res_def.jl:71
  [9] (::SlimOptim.var"#objgrad!#11"{var"#3#4", result{Float32}})(g::Matrix{Float32}, x::Matrix{Float32})
    @ SlimOptim ~/.julia/packages/SlimOptim/bTmY2/src/SPGSlim.jl:99
 [10] _spg(obj::Function, grad!::SlimOptim.var"#grad!#10"{var"#3#4", result{Float32}}, objgrad!::SlimOptim.var"#objgrad!#11"{var"#3#4", result{Float32}}, projection::SlimOptim.var"#projection#9"{typeof(proj), result{Float32}}, x::Matrix{Float32}, g::Matrix{Float32}, sol::result{Float32}, ls::Nothing, options::SlimOptim.SPG_params; callback::typeof(SlimOptim.noop_callback))
    @ SlimOptim ~/.julia/packages/SlimOptim/bTmY2/src/SPGSlim.jl:172
 [11] spg(funObj::var"#3#4", x::Matrix{Float32}, funProj::typeof(proj), options::SlimOptim.SPG_params; ls::Nothing, callback::Function)
    @ SlimOptim ~/.julia/packages/SlimOptim/bTmY2/src/SPGSlim.jl:103
 [12] spg(funObj::Function, x::Matrix{Float32}, funProj::Function, options::SlimOptim.SPG_params)
    @ SlimOptim ~/.julia/packages/SlimOptim/bTmY2/src/SPGSlim.jl:89
 [13] top-level scope
    @ /media/swina/WORK_DISK_1/Julia/FWI-Julia/scripts/02_res_def.jl:80
 [14] include(fname::String)
    @ Base.MainInclude ./client.jl:476
 [15] top-level scope
    @ REPL[2]:1
in expression starting at /media/swina/WORK_DISK_1/Julia/FWI-Julia/scripts/02_res_def.jl:80

Can anyone help me? I don't know how to solve these error!
Thanks

[mloubout]

What version of JUDI are you using? This is a new feature so you'll need the latest version.

[mloubout]

You will need JUDI 3.1.12 at the minimum

[asidosaputra]

Thank you, sir; I have updated my JUDI to the latest version, and working.