search

slimgroup / JUDI.jl

Julia Devito inversion.
https://slimgroup.github.io/JUDI.jl
MIT License
98 stars 31 forks source link

`DataPreconditioner` works strange #213

Closed kerim371 closed 6 months ago

kerim371 commented 11 months ago

Hi,

I'm trying to compute FWI using turning waves and filter source wavelet and data.

I have something similar to that:

# Data muters (t0 somehow can't be 0)
Ml_tur = judiDataMute(q.geometry, d_obs.geometry, vp=1500f0, t0=0.001f0, mode=:turning)    # keep turning waves

# Bandpass filter
Ml_freq = judiFilter(d_obs.geometry, 0.001, 3f0)
Mr_freq = judiFilter(q.geometry, 0.001, 3f0)

fval, gradient = fwi_objective(model0, Mr_freq[indsrc]*q[indsrc], Ml_tur[indsrc]*Ml_freq[indsrc]*get_data(d_obs[indsrc]), options=jopt, misfit=myloss)

I compared the result with:

fval, gradient = fwi_objective(model0, Mr_freq[indsrc]*q[indsrc], Ml_freq[indsrc]*get_data(d_obs[indsrc]), options=jopt, misfit=myloss)

and with:

fval, gradient = fwi_objective(model0, q[indsrc], get_data(d_obs[indsrc]), options=jopt, misfit=myloss)

Calculated gradient in these three cases stays the same. It looks like DataPreconditioners don't work in all of these three situations.

It only changes when using linear muter without bandpass filtering:

fval, gradient = fwi_objective(model0, q[indsrc], Ml_tur[indsrc]*get_data(d_obs[indsrc]), options=jopt, misfit=myloss)

Maybe I don't understand something?

As a conclusion:

  1. filtration doesn't work
  2. muting+filtration doesn't work
  3. muting alone works

JUDI v3.3.8

mloubout commented 11 months ago

What does the muted/filtered data look like compared to the clean one?

kerim371 commented 11 months ago

@mloubout here are the gradients (picture titles like 5.0 Hz doesn't mean anything): image

kerim371 commented 11 months ago

@mloubout I'm sorry, in my previous post I've sent slightly incorrect images.

I rechecked the problem and here is the update of information:

  1. muter doesn't work with filter (filter works and muter not)
  2. muter is not sensistive to what velocity I set: 1ms, 1500ms, 1e6m/s give the same gradient

image

mloubout commented 11 months ago

Before looking at those gradient

  1. I would first check what those preconditioners do on d_obs by themselves
  2. If you use those preconditioners with the lsrtm/fwi objectives you'll have to adapt the loss function so they get applied to the synthetic data as well. THis is after making sure that 1. does what you expect to the data
kerim371 commented 11 months ago

@mloubout I recheked the preconditionners info from preconditionners notebook example and some some understanding has come.

For example judiDataMute may work with negative t0 (#170) if we explicitly pass taperwidth:

Dmr = judiDataMute(q, d_obs, vp=1500f0, t0=-1f0; mode=:reflection, taperwidth=2)
Dmt = judiDataMute(q, d_obs, vp=1500f0, t0=-1f0; mode=:turning, taperwidth=2)

Also we have to explicitly pass taperwidth with very small t0 like 0.001 or there may be no muting. This is explained how taperwidth is defined within this function: taperwidth=floor(Int, 2/t0).

image

2. If you use those preconditioners with the lsrtm/fwi objectives you'll have to adapt the loss function so they get applied to the synthetic data as well. THis is after making sure that 1. does what you expect to the data

I guess I have to only implement muting there as filtration is done using source wavelet filtering. The muting probably should be done in the way similar to that:

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

but the problem is that Ml_tur is defined for all sources and typeof(dsyn) is a Matrix. Thus I can't get any information how to mute dsyn. The error is:

ERROR: TaskFailedException

    nested task error: JOLI.joUtilsException("*(jo,mvec) not implemented or type mismatch")
    Stacktrace:
     [1] jo_method_error(L::DataMute{Float32, :turning}, R::Matrix{Float32}, s::String)
       @ JOLI ~/Documents/Colada/r/julia-1.6/.julia/packages/JOLI/S3LJy/src/joUtils.jl:53
     [2] *(A::DataMute{Float32, :turning}, mv::Matrix{Float32})
       @ JOLI ~/Documents/Colada/r/julia-1.6/.julia/packages/JOLI/S3LJy/src/joAbstractOperator/base_functions.jl:74
     [3] myloss(dsyn::Matrix{Float32}, dobs::Matrix{Float32})
       @ Main /mnt/HDD_2TB/DMNG/model_1/fwi/fwi_spg.jl:180
     [4] top-level scope
       @ none:1
     [5] eval
       @ ./boot.jl:360 [inlined]
     [6] (::Distributed.var"#160#162"{Module, Expr})()
       @ Distributed ./task.jl:417
Stacktrace:
 [1] sync_end(c::Channel{Any})
   @ Base ./task.jl:369
 [2] macro expansion
   @ ./task.jl:388 [inlined]
 [3] remotecall_eval(m::Module, procs::Vector{Int64}, ex::Expr)
   @ Distributed /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/Distributed/src/macros.jl:223
 [4] top-level scope
   @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/Distributed/src/macros.jl:207
mloubout commented 11 months ago

dsyn will be sampled at get_dt(model) computational timestep rather than the "4ms" so you needa preconditionner setup for that time axis

kerim371 commented 11 months ago

@mloubout to set this up do I need to create a geometry with dt == get_dt(model) and pass these geometries to judiDataMute? I'm trying to do that in this way:

model_dt = get_dt(model0)
ntComp = get_computational_nt(q.geometry,d_obs.geometry,model0)
t_syn = model_dt*(ntComp[1]-1)
src_geom_syn = Geometry(container; key = "source", segy_depth_key = segy_depth_key_src, dt = model_dt, t = t_syn)
rec_geom_syn = Geometry(container; key = "receiver", segy_depth_key = segy_depth_key_rec, dt = model_dt, t = t_syn)

Ml_tur_syn = judiDataMute(src_geom_syn, rec_geom_syn, vp=1500f0, t0=0f0, mode=:turning, taperwidth=2)

function myloss(d_syn, d_obs)
  d_syn = Ml_tur_syn*d_syn
  fval = .5f0 * norm(d_syn - d_obs)^2
  adjoint_source = d_syn - d_obs
  return fval, adjoint_source
end

and I get an exception:

ERROR: TaskFailedException

    nested task error: JOLI.joUtilsException("*(jo,mvec) not implemented or type mismatch")
    Stacktrace:
     [1] jo_method_error(L::DataMute{Float32, :turning}, R::Matrix{Float32}, s::String)
       @ JOLI ~/Documents/Colada/r/julia-1.6/.julia/packages/JOLI/S3LJy/src/joUtils.jl:53
     [2] *(A::DataMute{Float32, :turning}, mv::Matrix{Float32})
       @ JOLI ~/Documents/Colada/r/julia-1.6/.julia/packages/JOLI/S3LJy/src/joAbstractOperator/base_functions.jl:74
     [3] myloss(d_syn::Matrix{Float32}, d_obs::Matrix{Float32})
       @ Main /mnt/HDD_2TB/DMNG/model_1/fwi/fwi_spg.jl:190
     [4] top-level scope
       @ none:1
     [5] eval
       @ ./boot.jl:360 [inlined]
     [6] (::Distributed.var"#160#162"{Module, Expr})()
       @ Distributed ./task.jl:417
Stacktrace:
 [1] sync_end(c::Channel{Any})
   @ Base ./task.jl:369
 [2] macro expansion
   @ ./task.jl:388 [inlined]
 [3] remotecall_eval(m::Module, procs::Vector{Int64}, ex::Expr)
   @ Distributed /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/Distributed/src/macros.jl:223
 [4] top-level scope
   @ /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/Distributed/src/macros.jl:207
mloubout commented 11 months ago

DataMute expect either a JUDI type (judiVector) or a pure vector as a standard mat-vec product so you need Ml_tur_syn*d_syn[:]

kerim371 commented 11 months ago

@mloubout I have changed this line as you recommended (d_syn = Ml_tur_syn*d_syn[:]) and the error stack now:

ERROR: (in a Julia function called from Python)
JULIA: BoundsError: attempt to access 4034×501×1 Array{Float32, 3} at index [1:4034, 1:501, 2]
Stacktrace:
  [1] throw_boundserror(A::Array{Float32, 3}, I::Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, Int64})
    @ Base ./abstractarray.jl:651
  [2] checkbounds
    @ ./abstractarray.jl:616 [inlined]
  [3] view(::Array{Float32, 3}, ::Function, ::Function, ::Int64)
    @ Base ./subarray.jl:177
  [4] muteshot(shot::Vector{Float32}, srcGeom::GeometryOOC{Float32}, recGeom::GeometryOOC{Float32}; vp::Vector{Float32}, t0::Vector{Float32}, mode::Symbol, taperwidth::Vector{Int64})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Preconditioners/DataPreconditioners.jl:111
  [5] matvec(D::DataMute{Float32, :turning}, x::Vector{Float32})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Preconditioners/DataPreconditioners.jl:64
  [6] *(J::DataMute{Float32, :turning}, v::Vector{Float32})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Preconditioners/base.jl:30
  [7] myloss(d_syn::Matrix{Float32}, d_obs::Matrix{Float32})
    @ Main /mnt/HDD_2TB/DMNG/model_1/fwi/fwi_spg.jl:191
  [8] (::PyCall.FuncWrapper{Tuple{Matrix{Float32}, Matrix{Float32}}, typeof(myloss)})(::Matrix{Float32}, ::Vararg{Matrix{Float32}, N} where N; kws::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/callback.jl:56
  [9] (::PyCall.FuncWrapper{Tuple{Matrix{Float32}, Matrix{Float32}}, typeof(myloss)})(::Matrix{Float32}, ::Vararg{Matrix{Float32}, N} where N)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/callback.jl:56
 [10] invokelatest(::Any, ::Any, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Base ./essentials.jl:708
 [11] invokelatest(::Any, ::Any, ::Vararg{Any, N} where N)
    @ Base ./essentials.jl:706
 [12] _pyjlwrap_call(f::PyCall.FuncWrapper{Tuple{Matrix{Float32}, Matrix{Float32}}, typeof(myloss)}, args_::Ptr{PyCall.PyObject_struct}, kw_::Ptr{PyCall.PyObject_struct})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/callback.jl:28
 [13] pyjlwrap_call(self_::Ptr{PyCall.PyObject_struct}, args_::Ptr{PyCall.PyObject_struct}, kw_::Ptr{PyCall.PyObject_struct})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/callback.jl:44
 [14] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:108 [inlined]
 [15] #107
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:43 [inlined]
 [16] disable_sigint
    @ ./c.jl:458 [inlined]
 [17] __pycall!
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:42 [inlined]
 [18] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, nargs::Int64, kw::PyCall.PyObject)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:29
 [19] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:11
 [20] #pycall#112
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:80 [inlined]
 [21] #4
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:82 [inlined]
 [22] (::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}}})()
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:74
 [23] lock(f::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}}}, l::ReentrantLock)
    @ Base ./lock.jl:187
 [24] pylock(f::JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:71
 [25] rlock_pycall(::PyCall.PyObject, ::Type{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}}, ::PyCall.PyObject, ::Vararg{Any, N} where N; kw::Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:81
 [26] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:50 [inlined]
 [27] macro expansion
    @ ./timing.jl:287 [inlined]
 [28] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:141 [inlined]
 [29] multi_src_fg(model_full::JUDI.IsoModel{Float32, 2}, source::judiVector{Float32, Matrix{Float32}}, dObs::judiVector{Float32, Matrix{Float32}}, dm::Nothing, options::JUDIOptions, nlind::Bool, lin::Bool, misfit::Function)
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:49
 [30] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:36 [inlined]
 [31] macro expansion
    @ ./timing.jl:287 [inlined]
 [32] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:141 [inlined]
 [33] run_and_reduce(func::Function, #unused#::Nothing, nsrc::Int64, arg_func::JUDI.var"#270#271"{JUDIOptions, Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:nlind, :lin, :misfit), Tuple{Bool, Bool, typeof(myloss)}}}, JUDI.IsoModel{Float32, 2}, judiVector{Float32, Matrix{Float32}}, judiVector{Float32, Matrix{Float32}}, Nothing})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:35
 [34] multi_src_fg!(G::PhysicalParameter{Float32, 2}, model::JUDI.IsoModel{Float32, 2}, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}, dm::Nothing; options::JUDIOptions, kw::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:nlind, :lin, :misfit), Tuple{Bool, Bool, typeof(myloss)}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:108
 [35] #multi_exp_fg!#249
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:191 [inlined]
 [36] fwi_objective!(G::PhysicalParameter{Float32, 2}, model::JUDI.IsoModel{Float32, 2}, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}; options::JUDIOptions, kw::Base.Iterators.Pairs{Symbol, typeof(myloss), Tuple{Symbol}, NamedTuple{(:misfit,), Tuple{typeof(myloss)}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:172
 [37] fwi_objective(model::JUDI.IsoModel{Float32, 2}, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}; options::JUDIOptions, kw::Base.Iterators.Pairs{Symbol, typeof(myloss), Tuple{Symbol}, NamedTuple{(:misfit,), Tuple{typeof(myloss)}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:137
 [38] objective_function(m_update::Matrix{Float64})
    @ Main /mnt/HDD_2TB/DMNG/model_1/fwi/fwi_spg.jl:225
 [39] (::SlimOptim.var"#objgrad!#11"{typeof(objective_function), result{Float32}})(g::Matrix{Float32}, x::Matrix{Float64})
    @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:99
 [40] _spg(obj::Function, grad!::SlimOptim.var"#grad!#10"{typeof(objective_function), result{Float32}}, objgrad!::SlimOptim.var"#objgrad!#11"{typeof(objective_function), 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 ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:172
 [41] spg(funObj::typeof(objective_function), x::Matrix{Float32}, funProj::typeof(proj), options::SlimOptim.SPG_params; ls::Nothing, callback::Function)
    @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:103
 [42] spg(funObj::Function, x::Matrix{Float32}, funProj::Function, options::SlimOptim.SPG_params)
    @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:92
 [43] top-level scope
    @ /mnt/HDD_2TB/DMNG/model_1/fwi/fwi_spg.jl:292
 [44] eval
    @ ./boot.jl:360 [inlined]
 [45] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1116
 [46] include_string(m::Module, txt::String, fname::String)
    @ Base ./loading.jl:1126
 [47] invokelatest(::Any, ::Any, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Base ./essentials.jl:708
 [48] invokelatest(::Any, ::Any, ::Vararg{Any, N} where N)
    @ Base ./essentials.jl:706
 [49] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/eval.jl:263
 [50] (::VSCodeServer.var"#65#70"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/eval.jl:181
 [51] withpath(f::VSCodeServer.var"#65#70"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/repl.jl:274
 [52] (::VSCodeServer.var"#64#69"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/eval.jl:179
 [53] hideprompt(f::VSCodeServer.var"#64#69"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/repl.jl:38
 [54] (::VSCodeServer.var"#63#68"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/eval.jl:150
 [55] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging ./logging.jl:491
 [56] with_logger
    @ ./logging.jl:603 [inlined]
 [57] (::VSCodeServer.var"#62#67"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/eval.jl:255
 [58] #invokelatest#2
    @ ./essentials.jl:708 [inlined]
 [59] invokelatest(::Any)
    @ Base ./essentials.jl:706
 [60] macro expansion
    @ ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/eval.jl:34 [inlined]
Stacktrace:
  [1] pyerr_check
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:75 [inlined]
  [2] pyerr_check
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:79 [inlined]
  [3] _handle_error(msg::String)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:96
  [4] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:110 [inlined]
  [5] #107
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:43 [inlined]
  [6] disable_sigint
    @ ./c.jl:458 [inlined]
  [7] __pycall!
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:42 [inlined]
  [8] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, nargs::Int64, kw::PyCall.PyObject)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:29
  [9] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:11
 [10] #pycall#112
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:80 [inlined]
 [11] #4
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:82 [inlined]
 [12] (::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}}})()
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:74
 [13] lock(f::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}}}, l::ReentrantLock)
    @ Base ./lock.jl:187
 [14] pylock(f::JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:71
 [15] rlock_pycall(::PyCall.PyObject, ::Type{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}}, ::PyCall.PyObject, ::Vararg{Any, N} where N; kw::Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:81
 [16] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:50 [inlined]
 [17] macro expansion
    @ ./timing.jl:287 [inlined]
 [18] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:141 [inlined]
 [19] multi_src_fg(model_full::JUDI.IsoModel{Float32, 2}, source::judiVector{Float32, Matrix{Float32}}, dObs::judiVector{Float32, Matrix{Float32}}, dm::Nothing, options::JUDIOptions, nlind::Bool, lin::Bool, misfit::Function)
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:49
 [20] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:36 [inlined]
 [21] macro expansion
    @ ./timing.jl:287 [inlined]
 [22] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:141 [inlined]
 [23] run_and_reduce(func::Function, #unused#::Nothing, nsrc::Int64, arg_func::JUDI.var"#270#271"{JUDIOptions, Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:nlind, :lin, :misfit), Tuple{Bool, Bool, typeof(myloss)}}}, JUDI.IsoModel{Float32, 2}, judiVector{Float32, Matrix{Float32}}, judiVector{Float32, Matrix{Float32}}, Nothing})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:35
 [24] multi_src_fg!(G::PhysicalParameter{Float32, 2}, model::JUDI.IsoModel{Float32, 2}, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}, dm::Nothing; options::JUDIOptions, kw::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:nlind, :lin, :misfit), Tuple{Bool, Bool, typeof(myloss)}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:108
 [25] #multi_exp_fg!#249
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:191 [inlined]
 [26] fwi_objective!(G::PhysicalParameter{Float32, 2}, model::JUDI.IsoModel{Float32, 2}, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}; options::JUDIOptions, kw::Base.Iterators.Pairs{Symbol, typeof(myloss), Tuple{Symbol}, NamedTuple{(:misfit,), Tuple{typeof(myloss)}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:172
 [27] fwi_objective(model::JUDI.IsoModel{Float32, 2}, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}; options::JUDIOptions, kw::Base.Iterators.Pairs{Symbol, typeof(myloss), Tuple{Symbol}, NamedTuple{(:misfit,), Tuple{typeof(myloss)}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:137
 [28] objective_function(m_update::Matrix{Float64})
    @ Main /mnt/HDD_2TB/DMNG/model_1/fwi/fwi_spg.jl:225
 [29] (::SlimOptim.var"#objgrad!#11"{typeof(objective_function), result{Float32}})(g::Matrix{Float32}, x::Matrix{Float64})
    @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:99
 [30] _spg(obj::Function, grad!::SlimOptim.var"#grad!#10"{typeof(objective_function), result{Float32}}, objgrad!::SlimOptim.var"#objgrad!#11"{typeof(objective_function), 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 ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:172
 [31] spg(funObj::typeof(objective_function), x::Matrix{Float32}, funProj::typeof(proj), options::SlimOptim.SPG_params; ls::Nothing, callback::Function)
    @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:103
 [32] spg(funObj::Function, x::Matrix{Float32}, funProj::Function, options::SlimOptim.SPG_params)
    @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:92
 [33] top-level scope
    @ /mnt/HDD_2TB/DMNG/model_1/fwi/fwi_spg.jl:292
mloubout commented 11 months ago

Looks liek you are trying to use the multi-source preconditioned on a single source shot record

kerim371 commented 11 months ago

@mloubout actually I have a 2D seismic line wich consists of a multiple shots. Thus geometry contains all these shots. If I understand you correctly I should have something like:

indsrc = 1

function myloss(d_syn, d_obs)
  d_syn = Ml_tur_syn[indsrc]*d_syn[:] # process only single shot [indsrc]
  fval = .5f0 * norm(d_syn - d_obs)^2
  adjoint_source = d_syn - d_obs
  return fval, adjoint_source
end

# process only single shot [indsrc]
fval, gradient = fwi_objective(model0, Mr_freq[indsrc]*q[indsrc], Ml_tur[indsrc]*Ml_freq[indsrc]*d_obs[indsrc], options=jopt, misfit=myloss)

In this case the error (it seems something should be reshaped):

ERROR: (in a Julia function called from Python)
JULIA: DimensionMismatch("dimensions must match: a has dims (Base.OneTo(4034), Base.OneTo(501)), b has dims (Base.OneTo(2021034),), mismatch at 1")
Stacktrace:
  [1] promote_shape
    @ ./indices.jl:178 [inlined]
  [2] promote_shape
    @ ./indices.jl:174 [inlined]
  [3] promote_shape(a::Vector{Float32}, b::Matrix{Float32})
    @ Base ./indices.jl:169
  [4] -(A::Vector{Float32}, B::Matrix{Float32})
    @ Base ./arraymath.jl:38
  [5] myloss(d_syn::Matrix{Float32}, d_obs::Matrix{Float32})
    @ Main /mnt/HDD_2TB/DMNG/model_1/fwi/fwi_spg.jl:192
  [6] (::PyCall.FuncWrapper{Tuple{Matrix{Float32}, Matrix{Float32}}, typeof(myloss)})(::Matrix{Float32}, ::Vararg{Matrix{Float32}, N} where N; kws::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/callback.jl:56
  [7] (::PyCall.FuncWrapper{Tuple{Matrix{Float32}, Matrix{Float32}}, typeof(myloss)})(::Matrix{Float32}, ::Vararg{Matrix{Float32}, N} where N)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/callback.jl:56
  [8] invokelatest(::Any, ::Any, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Base ./essentials.jl:708
  [9] invokelatest(::Any, ::Any, ::Vararg{Any, N} where N)
    @ Base ./essentials.jl:706
 [10] _pyjlwrap_call(f::PyCall.FuncWrapper{Tuple{Matrix{Float32}, Matrix{Float32}}, typeof(myloss)}, args_::Ptr{PyCall.PyObject_struct}, kw_::Ptr{PyCall.PyObject_struct})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/callback.jl:28
 [11] pyjlwrap_call(self_::Ptr{PyCall.PyObject_struct}, args_::Ptr{PyCall.PyObject_struct}, kw_::Ptr{PyCall.PyObject_struct})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/callback.jl:44
 [12] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:108 [inlined]
 [13] #107
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:43 [inlined]
 [14] disable_sigint
    @ ./c.jl:458 [inlined]
 [15] __pycall!
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:42 [inlined]
 [16] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, nargs::Int64, kw::PyCall.PyObject)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:29
 [17] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:11
 [18] #pycall#112
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:80 [inlined]
 [19] #4
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:82 [inlined]
 [20] (::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}}})()
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:74
 [21] lock(f::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}}}, l::ReentrantLock)
    @ Base ./lock.jl:187
 [22] pylock(f::JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:71
 [23] rlock_pycall(::PyCall.PyObject, ::Type{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}}, ::PyCall.PyObject, ::Vararg{Any, N} where N; kw::Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:81
 [24] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:50 [inlined]
 [25] macro expansion
    @ ./timing.jl:287 [inlined]
 [26] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:141 [inlined]
 [27] multi_src_fg(model_full::JUDI.IsoModel{Float32, 2}, source::judiVector{Float32, Matrix{Float32}}, dObs::judiVector{Float32, Matrix{Float32}}, dm::Nothing, options::JUDIOptions, nlind::Bool, lin::Bool, misfit::Function)
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:49
 [28] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:36 [inlined]
 [29] macro expansion
    @ ./timing.jl:287 [inlined]
 [30] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:141 [inlined]
 [31] run_and_reduce(func::Function, #unused#::Nothing, nsrc::Int64, arg_func::JUDI.var"#270#271"{JUDIOptions, Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:nlind, :lin, :misfit), Tuple{Bool, Bool, typeof(myloss)}}}, JUDI.IsoModel{Float32, 2}, judiVector{Float32, Matrix{Float32}}, judiVector{Float32, Matrix{Float32}}, Nothing})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:35
 [32] multi_src_fg!(G::PhysicalParameter{Float32, 2}, model::JUDI.IsoModel{Float32, 2}, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}, dm::Nothing; options::JUDIOptions, kw::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:nlind, :lin, :misfit), Tuple{Bool, Bool, typeof(myloss)}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:108
 [33] #multi_exp_fg!#249
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:191 [inlined]
 [34] fwi_objective!(G::PhysicalParameter{Float32, 2}, model::JUDI.IsoModel{Float32, 2}, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}; options::JUDIOptions, kw::Base.Iterators.Pairs{Symbol, typeof(myloss), Tuple{Symbol}, NamedTuple{(:misfit,), Tuple{typeof(myloss)}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:172
 [35] fwi_objective(model::JUDI.IsoModel{Float32, 2}, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}; options::JUDIOptions, kw::Base.Iterators.Pairs{Symbol, typeof(myloss), Tuple{Symbol}, NamedTuple{(:misfit,), Tuple{typeof(myloss)}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:137
 [36] objective_function(m_update::Matrix{Float64})
    @ Main /mnt/HDD_2TB/DMNG/model_1/fwi/fwi_spg.jl:225
 [37] (::SlimOptim.var"#objgrad!#11"{typeof(objective_function), result{Float32}})(g::Matrix{Float32}, x::Matrix{Float64})
    @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:99
 [38] _spg(obj::Function, grad!::SlimOptim.var"#grad!#10"{typeof(objective_function), result{Float32}}, objgrad!::SlimOptim.var"#objgrad!#11"{typeof(objective_function), 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 ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:172
 [39] spg(funObj::typeof(objective_function), x::Matrix{Float32}, funProj::typeof(proj), options::SlimOptim.SPG_params; ls::Nothing, callback::Function)
    @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:103
 [40] spg(funObj::Function, x::Matrix{Float32}, funProj::Function, options::SlimOptim.SPG_params)
    @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:92
 [41] top-level scope
    @ /mnt/HDD_2TB/DMNG/model_1/fwi/fwi_spg.jl:292
 [42] eval
    @ ./boot.jl:360 [inlined]
 [43] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1116
 [44] include_string(m::Module, txt::String, fname::String)
    @ Base ./loading.jl:1126
 [45] invokelatest(::Any, ::Any, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Base ./essentials.jl:708
 [46] invokelatest(::Any, ::Any, ::Vararg{Any, N} where N)
    @ Base ./essentials.jl:706
 [47] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/eval.jl:263
 [48] (::VSCodeServer.var"#65#70"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/eval.jl:181
 [49] withpath(f::VSCodeServer.var"#65#70"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/repl.jl:274
 [50] (::VSCodeServer.var"#64#69"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/eval.jl:179
 [51] hideprompt(f::VSCodeServer.var"#64#69"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/repl.jl:38
 [52] (::VSCodeServer.var"#63#68"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/eval.jl:150
 [53] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging ./logging.jl:491
 [54] with_logger
    @ ./logging.jl:603 [inlined]
 [55] (::VSCodeServer.var"#62#67"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/eval.jl:255
 [56] #invokelatest#2
    @ ./essentials.jl:708 [inlined]
 [57] invokelatest(::Any)
    @ Base ./essentials.jl:706
 [58] macro expansion
    @ ~/.vscode/extensions/julialang.language-julia-1.56.2/scripts/packages/VSCodeServer/src/eval.jl:34 [inlined]
Stacktrace:
  [1] pyerr_check
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:75 [inlined]
  [2] pyerr_check
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:79 [inlined]
  [3] _handle_error(msg::String)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:96
  [4] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/exception.jl:110 [inlined]
  [5] #107
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:43 [inlined]
  [6] disable_sigint
    @ ./c.jl:458 [inlined]
  [7] __pycall!
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:42 [inlined]
  [8] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, nargs::Int64, kw::PyCall.PyObject)
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:29
  [9] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}, kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}})
    @ PyCall ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:11
 [10] #pycall#112
    @ ~/Documents/Colada/r/julia-1.6/.julia/packages/PyCall/ilqDX/src/pyfncall.jl:80 [inlined]
 [11] #4
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:82 [inlined]
 [12] (::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}}})()
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:74
 [13] lock(f::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}}}, l::ReentrantLock)
    @ Base ./lock.jl:187
 [14] pylock(f::JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}}, PyCall.PyObject, Tuple{PyCall.PyObject, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}, Matrix{Float32}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:71
 [15] rlock_pycall(::PyCall.PyObject, ::Type{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}}, ::PyCall.PyObject, ::Vararg{Any, N} where N; kw::Base.Iterators.Pairs{Symbol, Any, NTuple{13, Symbol}, NamedTuple{(:t_sub, :space_order, :checkpointing, :freq_list, :ic, :is_residual, :born_fwd, :nlind, :dft_sub, :f0, :return_obj, :misfit, :illum), Tuple{Int64, Int64, Bool, Nothing, String, Bool, Bool, Bool, Int64, Float32, Bool, PyCall.PyObject, Bool}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:81
 [16] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:50 [inlined]
 [17] macro expansion
    @ ./timing.jl:287 [inlined]
 [18] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:141 [inlined]
 [19] multi_src_fg(model_full::JUDI.IsoModel{Float32, 2}, source::judiVector{Float32, Matrix{Float32}}, dObs::judiVector{Float32, Matrix{Float32}}, dm::Nothing, options::JUDIOptions, nlind::Bool, lin::Bool, misfit::Function)
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:49
 [20] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:36 [inlined]
 [21] macro expansion
    @ ./timing.jl:287 [inlined]
 [22] macro expansion
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/JUDI.jl:141 [inlined]
 [23] run_and_reduce(func::Function, #unused#::Nothing, nsrc::Int64, arg_func::JUDI.var"#270#271"{JUDIOptions, Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:nlind, :lin, :misfit), Tuple{Bool, Bool, typeof(myloss)}}}, JUDI.IsoModel{Float32, 2}, judiVector{Float32, Matrix{Float32}}, judiVector{Float32, Matrix{Float32}}, Nothing})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:35
 [24] multi_src_fg!(G::PhysicalParameter{Float32, 2}, model::JUDI.IsoModel{Float32, 2}, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}, dm::Nothing; options::JUDIOptions, kw::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:nlind, :lin, :misfit), Tuple{Bool, Bool, typeof(myloss)}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/propagation.jl:108
 [25] #multi_exp_fg!#249
    @ ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:191 [inlined]
 [26] fwi_objective!(G::PhysicalParameter{Float32, 2}, model::JUDI.IsoModel{Float32, 2}, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}; options::JUDIOptions, kw::Base.Iterators.Pairs{Symbol, typeof(myloss), Tuple{Symbol}, NamedTuple{(:misfit,), Tuple{typeof(myloss)}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:172
 [27] fwi_objective(model::JUDI.IsoModel{Float32, 2}, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}; options::JUDIOptions, kw::Base.Iterators.Pairs{Symbol, typeof(myloss), Tuple{Symbol}, NamedTuple{(:misfit,), Tuple{typeof(myloss)}}})
    @ JUDI ~/Documents/Colada/r/julia-1.6/.julia/dev/JUDI/src/TimeModeling/Modeling/misfit_fg.jl:137
 [28] objective_function(m_update::Matrix{Float64})
    @ Main /mnt/HDD_2TB/DMNG/model_1/fwi/fwi_spg.jl:225
 [29] (::SlimOptim.var"#objgrad!#11"{typeof(objective_function), result{Float32}})(g::Matrix{Float32}, x::Matrix{Float64})
    @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:99
 [30] _spg(obj::Function, grad!::SlimOptim.var"#grad!#10"{typeof(objective_function), result{Float32}}, objgrad!::SlimOptim.var"#objgrad!#11"{typeof(objective_function), 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 ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:172
 [31] spg(funObj::typeof(objective_function), x::Matrix{Float32}, funProj::typeof(proj), options::SlimOptim.SPG_params; ls::Nothing, callback::Function)
    @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:103
 [32] spg(funObj::Function, x::Matrix{Float32}, funProj::Function, options::SlimOptim.SPG_params)
    @ SlimOptim ~/Documents/Colada/r/julia-1.6/.julia/dev/SlimOptim/src/SPGSlim.jl:92
 [33] top-level scope
    @ /mnt/HDD_2TB/DMNG/model_1/fwi/fwi_spg.jl:292
mloubout commented 11 months ago 
indsrc = 1

function myloss(d_syn, d_obs)
  d_syn = reshape(Ml_tur_syn[indsrc]*d_syn[:], size(d_syn) # process only single shot [indsrc]
  fval = .5f0 * norm(d_syn - d_obs)^2
  adjoint_source = d_syn - d_obs
  return fval, adjoint_source
end
kerim371 commented 11 months ago

@mloubout right, thank you!

And is there any solution of sythetic muting for multisource FWI?

kerim371 commented 11 months ago

@mloubout I'm thinking about avoiding using fwi_objective and get data muted. According to the examples the gradient can be calculated using:

F = judiModeling(model0, q.geometry, d_obs.geometry; options=jopt)
J = judiJacobian(F(model0), q)
d0 = F(model0)*q[indsrc]
gradient = J' * Ml_tur * Ml_freq * (d0 - d_obs)

Then how to compute the residual?

mloubout commented 11 months ago

The residual is d0 - d_obs so that would work yes

kerim371 commented 11 months ago

@mloubout I tried to compute single shot using fwi_objective (as you proposed previously) and gradient = J' * Ml_tur * Ml_freq * (d0 - d_obs). There is difference in gradient that is perhaps caused by filtration: if I don't use frequency filtering then the gradient is equal. image

kerim371 commented 11 months ago

@mloubout I'm trying to understand the difference between using operators and fwi_objective. For example these two snippets give exactly same gradient:

jopt = JUDI.Options(
    subsampling_factor=10,
    IC = "fwi",
    limit_m = true,
    buffer_size = buffer_size,
    optimal_checkpointing=false,
    free_surface=false, 
    space_order=8)

F = judiModeling(model0, q.geometry, d_obs.geometry; options=jopt)
J = judiJacobian(F(model0), q)

indsrc = 1
d0 = F(model0)[indsrc]*q[indsrc]
d = d_obs[indsrc]
r = d0 - d
gradient = J'[indsrc] * r

and fval, gradient = fwi_objective(model0, q[indsrc], d_obs[indsrc], options=jopt). Both these snippets of code produce the same gradient.

But if I add frequency filtering the gradients becomes different like I showed in the post above:

Ml_freq = judiFilter(d_obs.geometry, 0.001, 2f0)
Mr_freq = judiFilter(q.geometry, 0.001, 2f0)

d0 = F(model0)[indsrc]*(Mr_freq[indsrc]*q[indsrc])
d = Ml_freq[indsrc]*d_obs[indsrc]
r = d0 - d
gradient = J'[indsrc] * r

and fval, gradient = fwi_objective(model0, Mr_freq[indsrc]*q[indsrc], Ml_freq[indsrc]*d_obs[indsrc], options=jopt).

I believe I'm missing something when doing frequency filtering on operators even if I follow the preconditionners example.

mloubout commented 11 months ago

You need to give the filtered source to J i.e J = judiJacobian(F(model0), Mr_freq*q) or the Jacobian will use the full band source for the forward wavefield

kerim371 commented 11 months ago

You need to give the filtered source to J i.e J = judiJacobian(F(model0), Mr_freq*q) or the Jacobian will use the full band source for the forward wavefield

Right! thank you very much!