Muting in FWI user provided misfit function.

[NikoYura]

Hello!

I want to add built-in JUDI muting (JUDI.DataMute) into user provided misfit function (an example):

using FFTW, Zygote

n2(x) = x / norm(x, 2)

function Hilbert(x)
    n = size(x, 1)
    σ = sign.(-n/2+1:n/2)
    y = imag(ifft(fftshift(σ.*fftshift(fft(x, 1), 1), 1), 1))
    return y
end

HLoss(dsyn, dobs) = sum(abs2.((dsyn - dobs) .+ 1im .* Hilbert(dsyn - dobs)))

function envelope(dsyn, dobs; Ml)
    #dsyn = Ml*dsyn
    ϕ = HLoss(n2(dsyn), n2(dobs))
    g = gradient(xs -> HLoss(n2(xs), n2(dobs)), dsyn)
    # Zygote always returns a tuple, in this case (g,)
    return ϕ, real.(g[1])
end

g_const = 0
opt = Options(limit_m = true, buffer_size = 3000f0, optimal_checkpointing = true)
function objective_function(x, d_obs, misfit=mse)
    model0.m .= reshape(x,model0.n);
    
    Ml = judiDataMute(q.geometry, d_obs.geometry)
    
    # fwi function value and gradient
    i = randperm(d_obs.nsrc)[1:batchsize]
    if misfit == envelope
        fval, grad = fwi_objective(model0, q[i], d_obs[i]; misfit=envelope(dsyn=missing, dobs=missing; Ml=Ml[i]), options = opt)
    else
        fval, grad = fwi_objective(model0, q[i], d_obs[i]; misfit=misfit, options = opt)
    end
    # Normalize for nicer convergence
    g = grad.data
    # Mute water
    #g[:, 1:15] .= 0f0
    g_const == 0 && (global g_const = 1/norm(g, Inf))
    return fval, g_const .* g
end

How to do it correctly because it doesn't work for me. Should I go into fwi_objective directly and make muting implicitly?

[mloubout]

Hi

Happy you're using JUDI.

So this is one of those corner case ones that are a bit tricky to get around and I don't have a clean solution at the moment and there is mainly two reasons

• the misfit is the same for all sources, so there would need to be a way to make it "per-source" so that it uses the per-source preconditioned with the right index
• The misift is computed during propagation, i.e on the computational time axis wich is get_dt(model) instead of the data sampling rate so there would be a need to be able to project the precondionner onto this time axis. This part might be not too complicated though it's just a dt change, as long as the previous one is done.

So the way to go would be to make a misfit "list"

function misifit(x, y, Ml)
   mldt = project(Ml, dt_scale=inferdt) # should be roughly size(x, 1)/Ml.geometry.nt[1]
   ...

• misifit = [(x, y)->misifit(x, y, Ml[i]) for i=1:nsrc]
• Add a case to get_src in src/TimeModeling/Modeling/propagation.jl to process it and scatter it to each task
• implement project(ml::DataMute, dt) that creates a new one with the new dt

Hope this helps

I am not sure when I'll have some time to spend on it though.