Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DataPreconditioner works strange #213

Closed
kerim371 opened this issue Nov 7, 2023 · 19 comments · Fixed by #232
Closed

DataPreconditioner works strange #213

kerim371 opened this issue Nov 7, 2023 · 19 comments · Fixed by #232

Comments

@kerim371
Copy link
Contributor

kerim371 commented Nov 7, 2023

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
Copy link
Member

mloubout commented Nov 7, 2023

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

@kerim371
Copy link
Contributor Author

kerim371 commented Nov 7, 2023

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

@kerim371
Copy link
Contributor Author

kerim371 commented Nov 7, 2023

@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
Copy link
Member

mloubout commented Nov 7, 2023

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
Copy link
Contributor Author

kerim371 commented Nov 7, 2023

@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
Copy link
Member

mloubout commented Nov 7, 2023

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

@kerim371
Copy link
Contributor Author

kerim371 commented Nov 7, 2023

@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
Copy link
Member

mloubout commented Nov 7, 2023

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
Copy link
Contributor Author

kerim371 commented Nov 7, 2023

@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
Copy link
Member

mloubout commented Nov 7, 2023

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

@kerim371
Copy link
Contributor Author

kerim371 commented Nov 7, 2023

@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

@kerim371 kerim371 closed this as completed Nov 7, 2023
@kerim371 kerim371 reopened this Nov 7, 2023
@mloubout
Copy link
Member

mloubout commented Nov 7, 2023

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
Copy link
Contributor Author

kerim371 commented Nov 7, 2023

@mloubout right, thank you!

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

@kerim371
Copy link
Contributor Author

kerim371 commented Nov 7, 2023

@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
Copy link
Member

mloubout commented Nov 7, 2023

The residual is d0 - d_obs so that would work yes

@kerim371
Copy link
Contributor Author

kerim371 commented Nov 7, 2023

@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
Copy link
Contributor Author

kerim371 commented Nov 8, 2023

@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
Copy link
Member

mloubout commented Nov 8, 2023

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
Copy link
Contributor Author

kerim371 commented Nov 8, 2023

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!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants