Trying to run FWI with pqn from SlimOptim

[kerim371]

Hi,

I started to test pqn algorithm and I rely on this tutorial

First of all there should be a typo: function defined as function objective(F0, G, m, dobs, wb_mask) but it is called with g(x) = objective(F0, x, dobs, wb_mask). As you can see G is missing. I guess we don't need it in function definition right?

But the problem is at first iteration objective(...) function brings m with somehow modified values: it has negative values (minimum(m) gives some negative number) but my starting model contains only positive model0.m.

It would be great if you could give some hints because for now there is a great number of optimization parameters that I don't know yet (like TV).

For the time economy I use only a single shot (that is why shot_from == shot_to) and the code looks like:

# SETUP CONSTARAINTS
options=PARSDMM_options()
options.FL=Float32
options=default_PARSDMM_options(options,options.FL)
constraint = Vector{SetIntersectionProjection.set_definitions}()
constraint2 = Vector{SetIntersectionProjection.set_definitions}()

# bounds:
set_type = "bounds"
TD_OP = "identity"
app_mode = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,vec(mminArr),vec(mmaxArr),app_mode,custom_TD_OP));
push!(constraint2, set_definitions(set_type,TD_OP,vec(mminArr),vec(mmaxArr),app_mode,custom_TD_OP));

#TV
(TV,dummy1,dummy2,dummy3) = get_TD_operator(model0,"TV",options.FL)
m_min = 0.0
m_max = norm(TV*vec(m0),1) * .5
set_type = "l1"
TD_OP = "TV"
app_mode = ("matrix","")
custom_TD_OP = ([],false)
push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP));

# set up constraints with bounds only, precompute some things and define projector
(P_sub2,TD_OP2,set_Prop2) = setup_constraints(constraint2, model0,options.FL)
(TD_OP2,AtA2,l2,y2) = PARSDMM_precompute_distribute(TD_OP2,set_Prop2,model0,options)
options2 = deepcopy(options)
options2.rho_ini = ones(length(TD_OP2))*10.0

proj_intersection2 = x-> PARSDMM(x, AtA2, TD_OP2, set_Prop2, P_sub2, model0, options2)  

# Projection function
function prj2(input)
    input = Float32.(input)
    (x,dummy1,dummy2,dymmy3) = proj_intersection2(vec(input))
    return x
end

# set up constraints with bounds and TV
(P_sub,TD_OP,set_Prop) = setup_constraints(constraint, model0,options.FL)
(TD_OP,AtA,l,y) = PARSDMM_precompute_distribute(TD_OP,set_Prop,model0,options)
options.rho_ini = ones(length(TD_OP))*10.0

proj_intersection = x-> PARSDMM(x, AtA, TD_OP, set_Prop, P_sub, model0, options)

# Projection function
function prj(input)
    input = Float32.(input)
    (x,dummy1,dummy2,dymmy3) = proj_intersection(vec(input))
    return x
end

############# FWI #############
# JUDI options
global jopt = JUDI.Options(
    IC = "fwi",
    limit_m = true,
    buffer_size = buffer_size,
    optimal_checkpointing=false,
    free_surface=false,  # free_surface is ON to model multiples as well
    space_order=8)     # increase space order for > 12 Hz source wavelet

F0 = judiModeling(model0, q.geometry, d_obs.geometry; options=jopt)

# build Jacoian
J = judiJacobian(F0, q)

wb_mask = ones(Float32,size(m0))
wb_mask[:, 1:seabed_ind] .= 0f0

# optimization parameters
niterations = 15
shot_from = 1
shot_to = 1
shot_step = 1   # we may want to calculate only each Nth shot to economy time
count = 0
fhistory = Vector{Float32}(undef, 0)
mute_reflections = false
mute_turning = false

function objective(m, d_obs, wb_mask)

  global x, z, count, jopt, seabed_ind;
  count += 1

  @info "minimum(m): $(minimum(m))"  # HERE THE MINIMAL VALUE IS NEGATIVE. I THINK AT FIRST ITERATION SHOULD BE model0.m == m ISN'T IT?
  model0.m .= Float32.(reshape(m, model0.n))
  t = @elapsed begin
    indsrc = shot_from:shot_step:shot_to
    if mute_reflections
        fval, gradient = fwi_objective(model0, Mr_freq[indsrc]*q[indsrc], Ml_tur[indsrc]*Ml_freq[indsrc]*get_data(d_obs[indsrc]), options=jopt)
    elseif mute_turning
        fval, gradient = fwi_objective(model0, Mr_freq[indsrc]*q[indsrc], Ml_ref[indsrc]*Ml_freq[indsrc]*get_data(d_obs[indsrc]), options=jopt)
    else
        fval, gradient = fwi_objective(model0, Mr_freq[indsrc]*q[indsrc], Ml_freq[indsrc]*get_data(d_obs[indsrc]), options=jopt)
    end
    @info "gradient: $gradient"
    gradient = reshape(gradient, model0.n)
  end
  gradient .*= wb_mask
  ϕ = .5*norm(d0 .- d_obs)^2
  if gscale == 0.0
      # compute scalar from first gradient, apply to future gradients
      global gscale = .25 ./ maximum(gradient) 
      @show gscale
  end
  gradient .*= gscale

  return ϕ, vec(gradient.data)
end

# struct to save the first gradient scalar
gscale = 0f0
g(x) = objective(x, d_obs, wb_mask)

# FWI with PQN
gscale = 0f0
options_pqn = pqn_options(progTol=0, store_trace=true, verbose=3, maxIter=niterations)
sol = pqn(g, vec(m0), prj, options_pqn);

[mloubout]

HERE THE MINIMAL VALUE IS NEGATIVE. I THINK AT FIR

This will be prj(vec(m0)) SO would double check what you have in m0 and mminArr/mmaxArr and make sure the prj does what you expect it to do first by itself.

[kerim371]

What actually prj should do?
Any information or link to what it is about may help

[mloubout]

It's a projection it imposes the constrains to the input. SO for example if you use bounds constraints then prj(x) return an array that has all its entries between the minimum and maximum bound. You can read the SetIntersectionProjection.jl package documentation for for info about the different projections.

[kerim371]

The problem is connected with push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP));
After I commented it the computations began but then I got the error:

ERROR: MethodError: no method matching linesearch(::BackTracking{Float64, Int64}, ::result{Float64}, ::Vector{Float32}, ::SlimOptim.var"#obj#24"{SlimOptim.var"#objgrad!#23"{typeof(g), result{Float64}}, Vector{Float64}}, ::SlimOptim.var"#grad!#22"{typeof(g), result{Float64}}, ::SlimOptim.var"#objgrad!#23"{typeof(g), result{Float64}}, ::Float64, ::Float64, ::Float64, ::Vector{Float64})
Closest candidates are:
  linesearch(::Any, ::result, ::AbstractArray{T, N} where N, ::Function, ::Function, ::Function, ::Any, ::Any, ::Any, ::AbstractArray{T, N} where N) where T at /home/kerim/Documents/Colada/r/julia-1.6/.julia/packages/SlimOptim/33Fm9/src/linesearches.jl:7
  linesearch(::Any, ::result, ::AbstractArray{T, N} where N, ::Function, ::Function, ::Function, ::T, ::T, ::T, ::AbstractArray{T, N} where N) where T at /home/kerim/Documents/Colada/r/julia-1.6/.julia/packages/SlimOptim/33Fm9/src/linesearches.jl:27
Stacktrace:
 [1] _pqn(obj::Function, grad!::SlimOptim.var"#grad!#22"{typeof(g), result{Float64}}, objgrad!::SlimOptim.var"#objgrad!#23"{typeof(g), result{Float64}}, projection::SlimOptim.var"#projection#21"{typeof(prj), result{Float64}}, x::Vector{Float64}, g::Vector{Float64}, sol::result{Float64}, ls::Nothing, options::SlimOptim.PQN_params; callback::typeof(SlimOptim.noop_callback))
   @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/packages/SlimOptim/33Fm9/src/PQNSlim.jl:262
 [2] pqn(funObj::typeof(g), x::Vector{Float64}, funProj::typeof(prj), options::SlimOptim.PQN_params; ls::Nothing, callback::Function)
   @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/packages/SlimOptim/33Fm9/src/PQNSlim.jl:95
 [3] pqn(funObj::Function, x::Vector{Float64}, funProj::Function, options::SlimOptim.PQN_params)
   @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/packages/SlimOptim/33Fm9/src/PQNSlim.jl:85
 [4] top-level scope
   @ /mnt/HDD_2TB/DMNG/model_1/fwi/fwi_pqn.jl:325

[mloubout]

You have bad matching types, for example you. use .25 instead of .25f0 for the gradient scaling so the datatypes won't be matching.

[kerim371]

I changed all the numbers to 1f0 but I still get printed info that BackTracking is Float64:

Line search: BackTracking{Float64, Int64}

Is there a way to set it to 32?

[mloubout]

It infers its type from x so your vec(m0), so your m0 is probably FLoat64

[kerim371]

Yes, completely right! thank you.
I will try now to solve other issues with boundary containts