search

slimgroup / JUDI.jl

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

Custom loss function without `fwi_objective` #230

Open kerim371 opened 4 months ago

kerim371 commented 4 months ago

Hi,

As far as I know fwi_objective doesn't support muting yet #213. Thus I do the FWI without fwi_objective. How to use custom loss function in that case? Or where in JUDI source code the operator - is defined for residual that I use for now (d0 - Ml_freq[indsrc]*d_obs[indsrc]) so I can override it?

@everywhere envelope(x, y) = sum(abs2.((x - y) .+ 1im .* H(x - y)))
@everywhere denvelope(x, y) = Zygote.gradient(xs->envelope(xs, y), x)[1]
@everywhere myloss(x, y) = (envelope(x, y), denvelope(x, y))

d0 = F(model0)[indsrc]*(Mr_freq[indsrc]*q[indsrc])
r = Ml_tur[indsrc] * (d0 - Ml_freq[indsrc]*d_obs[indsrc])
mloubout commented 4 months ago

JUDI objects (judiVector, judiWeights, judiWavefields) are all designed to behave as vectors so you can directly use the custom misfit on the judiVectors i.e

phi, r = myloss(Ml_tur[indsrc] * d0 , Ml_tur[indsrc] * d0 Ml_freq[indsrc]*d_obs[indsrc])

kerim371 commented 4 months ago

@mloubout it seems it doesn't work as is:

@everywhere envelope(x, y) = sum(abs2.((x - y) .+ 1im .* H(x - y)))
@everywhere denvelope(x, y) = Zygote.gradient(xs->envelope(xs, y), x)[1]
@everywhere myloss(x, y) = (envelope(x, y), denvelope(x, y))

d0 = F(model0)[indsrc]*(Mr_freq[indsrc]*q[indsrc])
phi, r = myloss(Ml_tur[indsrc] * d0 , Ml_tur[indsrc] * Ml_freq[indsrc]*d_obs[indsrc])

produces error:

ERROR: MethodError: Cannot `convert` an object of type judiVector{Float32, Matrix{Float32}} to an object of type ComplexF32

Closest candidates are:
  convert(::Type{T}, ::Base.TwicePrecision) where T<:Number
   @ Base twiceprecision.jl:273
  convert(::Type{T}, ::AbstractChar) where T<:Number
   @ Base char.jl:185
  convert(::Type{T}, ::CartesianIndex{1}) where T<:Number
   @ Base multidimensional.jl:127
  ...

Stacktrace:
  [1] setindex!(A::Vector{ComplexF32}, x::judiVector{Float32, Matrix{Float32}}, i1::Int64)
    @ Base ./array.jl:969
  [2] copyto_unaliased!(deststyle::IndexLinear, dest::Vector{ComplexF32}, srcstyle::IndexCartesian, src::judiVector{Float32, Matrix{Float32}})
    @ Base ./abstractarray.jl:1099
  [3] copyto!(dest::Vector{ComplexF32}, src::judiVector{Float32, Matrix{Float32}})
    @ Base ./abstractarray.jl:1073
  [4] circcopy!(dest::Vector{ComplexF32}, src::judiVector{Float32, Matrix{Float32}})
    @ Base ./multidimensional.jl:1244
  [5] copy1(#unused#::Type{ComplexF32}, x::judiVector{Float32, Matrix{Float32}})
    @ AbstractFFTs ~/Documents/Colada/r/julia-1.9/.julia/packages/AbstractFFTs/4iQz5/src/definitions.jl:54
  [6] complexfloat
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/AbstractFFTs/4iQz5/src/definitions.jl:47 [inlined]
  [7] fft
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/AbstractFFTs/4iQz5/src/definitions.jl:214 [inlined]
  [8] H(x::judiVector{Float32, Matrix{Float32}})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:184
  [9] envelope(x::judiVector{Float32, Matrix{Float32}}, y::judiVector{Float32, Matrix{Float32}})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:188
 [10] myloss(x::judiVector{Float32, Matrix{Float32}}, y::judiVector{Float32, Matrix{Float32}})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:192
 [11] objective_function(m_update::Matrix{Float64})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:238
 [12] (::SlimOptim.var"#objgrad!#11"{typeof(objective_function), result{Float32}})(g::Matrix{Float32}, x::Matrix{Float64})
    @ SlimOptim ~/Documents/Colada/r/julia-1.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:99
 [13] _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.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:172
 [14] 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.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:103
 [15] spg(funObj::Function, x::Matrix{Float32}, funProj::Function, options::SlimOptim.SPG_params)
    @ SlimOptim ~/Documents/Colada/r/julia-1.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:89
 [16] top-level scope
    @ ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:309

I tried to overload myloss function to take judiVector{Float32, Matrix{Float32}} as arguments (like this: @everywhere envelope(x::JUDI.judiVector{Float32, Matrix{Float32}}, y::JUDI.judiVector{Float32, Matrix{Float32}}) = sum([abs2.((xd - yd) .+ 1im * H(xd - yd)) for (xd,yd) in zip(x.data,y.data)])) but couldn't do the same with denvelope. I think I led astray at some point.

JUDI: v3.3.8

kerim371 commented 4 months ago

Could you please take a look.

This works fine, is it the optimal solution or there are other approaches?

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

@everywhere envelope(x, y) = sum(abs2.((x - y) .+ 1im .* H(x - y)))
@everywhere envelope(x::JUDI.judiVector{Float32, Matrix{Float32}}, y::JUDI.judiVector{Float32, Matrix{Float32}}) = sum([abs2.((xd - yd) .+ 1im * H(xd - yd)) for (xd,yd) in zip(x.data,y.data)])
@everywhere denvelope(x, y) = Zygote.gradient(xs->envelope(xs, y), x)[1]
@everywhere denvelope(x::JUDI.judiVector{Float32, Matrix{Float32}}, y::JUDI.judiVector{Float32, Matrix{Float32}}) = [Zygote.gradient(xs->envelope(xs, yd), xd)[1] for (xd,yd) in zip(x.data,y.data)]
@everywhere myloss(x, y) = (envelope(x, y), denvelope(x, y))

d0 = F(model0)[indsrc]*(Mr_freq[indsrc]*q[indsrc])
r = judiVector(d0.geometry, myloss(Ml_tur[indsrc] * d0 , Ml_tur[indsrc] * Ml_freq[indsrc]*d_obs[indsrc])[2])
gradient = J'[indsrc] * r
mloubout commented 4 months ago

No it needs the Complex convert to be defined this is just a local workaround it doesn't fix the issue itself

kerim371 commented 3 months ago

@mloubout hi,

I believe to perform fft(d_obs, 1) we need to add judi vector convertor to complex space:

Base.convert(::Type{Complex{T}}, v::JUDI.judiVector{T, Array{T, N}}) where {T<:Number, N} = judiVector(v.geometry, convert.(Complex{T}, v.data[1]))

But it seems that judi vector can't be constructed with complex array:

ERROR: MethodError: no method matching tof32(::Matrix{ComplexF32})

Closest candidates are:
  tof32(::Number)
   @ JUDI ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/JEsVr/src/TimeModeling/Types/abstract.jl:208
  tof32(::Array{Array{T, N}, 1}) where {N, T<:Real}
   @ JUDI ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/JEsVr/src/TimeModeling/Types/abstract.jl:210
  tof32(::StepRangeLen)
   @ JUDI ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/JEsVr/src/TimeModeling/Types/abstract.jl:212
  ...

Stacktrace:
  [1] judiVector(geometry::GeometryIC{Float32}, data::Matrix{ComplexF32})
    @ JUDI ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/JEsVr/src/TimeModeling/Types/judiVector.jl:77
  [2] convert(#unused#::Type{ComplexF32}, v::judiVector{Float32, Matrix{Float32}})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:201
  [3] setindex!(A::Vector{ComplexF32}, x::judiVector{Float32, Matrix{Float32}}, i1::Int64)
    @ Base ./array.jl:969
  [4] copyto_unaliased!(deststyle::IndexLinear, dest::Vector{ComplexF32}, srcstyle::IndexCartesian, src::judiVector{Float32, Matrix{Float32}})
    @ Base ./abstractarray.jl:1099
  [5] copyto!(dest::Vector{ComplexF32}, src::judiVector{Float32, Matrix{Float32}})
    @ Base ./abstractarray.jl:1073
  [6] circcopy!(dest::Vector{ComplexF32}, src::judiVector{Float32, Matrix{Float32}})
    @ Base ./multidimensional.jl:1244
  [7] copy1(#unused#::Type{ComplexF32}, x::judiVector{Float32, Matrix{Float32}})
    @ AbstractFFTs ~/Documents/Colada/r/julia-1.9/.julia/packages/AbstractFFTs/4iQz5/src/definitions.jl:54
  [8] complexfloat
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/AbstractFFTs/4iQz5/src/definitions.jl:47 [inlined]
  [9] fft
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/AbstractFFTs/4iQz5/src/definitions.jl:214 [inlined]
 [10] H(x::judiVector{Float32, Matrix{Float32}})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:208
 [11] envelope(x::judiVector{Float32, Matrix{Float32}}, y::judiVector{Float32, Matrix{Float32}})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:212
 [12] myloss(x::judiVector{Float32, Matrix{Float32}}, y::judiVector{Float32, Matrix{Float32}})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:216
 [13] objective_function(m_update::Matrix{Float64})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:262
 [14] (::SlimOptim.var"#objgrad!#11"{typeof(objective_function), result{Float32}})(g::Matrix{Float32}, x::Matrix{Float64})
    @ SlimOptim ~/Documents/Colada/r/julia-1.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:99
 [15] _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.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:172
 [16] 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.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:103
 [17] spg(funObj::Function, x::Matrix{Float32}, funProj::Function, options::SlimOptim.SPG_params)
    @ SlimOptim ~/Documents/Colada/r/julia-1.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:89
 [18] top-level scope
    @ ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:333

judiVector construcotr tries to convert complex numbers to Float32. Do we need to add support for complex judi vector?

P.S. sorry for such long reply

mloubout commented 3 months ago

You probably need to either

Base.convert(::Type{Complex{T}}, v::JUDI.judiVector{T, Array{T, N}}) where {T<:Number, N} = judiVector{Complex{T}, Matrix{Complex{T}}}(v.nsrc, v.geometry, convert.(Complex{T}, v.data[1]))

or add the Complex{T} case to tof32

kerim371 commented 3 months ago

This partly solves the problem:

Base.convert(::Type{Complex{T}}, v::JUDI.judiVector{T, Array{T, N}}) where {T<:Number, N} judiVector{Complex{T}, Matrix{Complex{T}}}(v.nsrc, v.geometry, [convert.(Complex{T}, v.data[1])])

but fft(d_obs, 1) gives the error:

ERROR: TypeError: in typeassert, expected ComplexF32, got a value of type judiVector{ComplexF32, Matrix{ComplexF32}}
Stacktrace:
 [1] setindex!(A::Vector{ComplexF32}, x::judiVector{Float32, Matrix{Float32}}, i1::Int64)
   @ Base ./array.jl:969
 [2] copyto_unaliased!(deststyle::IndexLinear, dest::Vector{ComplexF32}, srcstyle::IndexCartesian, src::judiVector{Float32, Matrix{Float32}})
   @ Base ./abstractarray.jl:1099
 [3] copyto!(dest::Vector{ComplexF32}, src::judiVector{Float32, Matrix{Float32}})
   @ Base ./abstractarray.jl:1073
 [4] circcopy!(dest::Vector{ComplexF32}, src::judiVector{Float32, Matrix{Float32}})
   @ Base ./multidimensional.jl:1244
 [5] copy1(#unused#::Type{ComplexF32}, x::judiVector{Float32, Matrix{Float32}})
   @ AbstractFFTs ~/Documents/Colada/r/julia-1.9/.julia/packages/AbstractFFTs/4iQz5/src/definitions.jl:54
 [6] complexfloat
   @ ~/Documents/Colada/r/julia-1.9/.julia/packages/AbstractFFTs/4iQz5/src/definitions.jl:47 [inlined]
 [7] fft(x::judiVector{Float32, Matrix{Float32}}, region::Int64)
   @ AbstractFFTs ~/Documents/Colada/r/julia-1.9/.julia/packages/AbstractFFTs/4iQz5/src/definitions.jl:214
 [8] top-level scope
   @ REPL[19]:1

seems to be I don't understand what exactly should return Base.convert(::Type{Complex{T}}, v::JUDI.judiVector{T, Array{T, N}})...

mloubout commented 3 months ago

Ok so thinking about it it's a bit tricky because of all the internals of FFTW so your best bet might be to just define

function fft(x::judiVector, dims)

  do ffft on data for dim
  return judiVector{Complex{T}, Matrix{Complex{T}}} ...
 end
kerim371 commented 3 months ago

Thank you! That works:

function JUDI.FFTW.fft(x::JUDI.judiVector{T, Array{T, N}}, dims) where {T<:Number, N}
    return judiVector{Complex{T}, Matrix{Complex{T}}}(x.nsrc, x.geometry, [fft(x.data[1], dims)])
end
mloubout commented 3 months ago

That will only work for single source if you use x.data[1] you need to do fft for all sources

kerim371 commented 3 months ago

Yes right, thank you!

function JUDI.FFTW.fft(x::JUDI.judiVector{T, Array{T, N}}, dims) where {T<:Number, N}
    return judiVector{Complex{T}, Matrix{Complex{T}}}(x.nsrc, x.geometry, [fft(v, dims) for v in x.data])
end

function JUDI.FFTW.ifft(x::JUDI.judiVector{T, Array{T, N}}, dims) where {T<:Number, N}
    return judiVector{T, Matrix{T}}(x.nsrc, x.geometry, [ifft(v, dims) for v in x.data])
end

@everywhere function H(x)
    n = size(x, 1)
    σ = Float32.(ifftshift(sign.(-n/2+1:n/2)))
    y = imag(ifft(σ.*fft(x, 1), 1))
    return y
end

@everywhere envelope(x, y) = sum(abs2.((x - y) .+ 1im .* H(x - y)))
@everywhere denvelope(x, y) = Zygote.gradient(xs->envelope(xs, y), x)[1]
@everywhere myloss(x, y) = (envelope(x, y), denvelope(x, y))

Now during the inversion I get the error:

ERROR: MethodError: no method matching length(::JUDI.MultiSource)

Closest candidates are:
  length(::Union{Base.KeySet, Base.ValueIterator})
   @ Base abstractdict.jl:58
  length(::Union{SparseArrays.FixedSparseVector{Tv, Ti}, SparseArrays.SparseVector{Tv, Ti}} where {Tv, Ti})
   @ SparseArrays ~/Documents/Colada/r/julia-1.9/share/julia/stdlib/v1.9/SparseArrays/src/sparsevector.jl:95
  length(::Union{JLD2.Group, JLD2.JLDFile})
   @ JLD2 ~/Documents/Colada/r/julia-1.9/.julia/packages/JLD2/oYEEg/src/JLD2.jl:500
  ...

Stacktrace:
  [1] _similar_shape(itr::JUDI.MultiSource, #unused#::Base.HasLength)
    @ Base ./array.jl:658
  [2] _collect(cont::UnitRange{Int64}, itr::JUDI.MultiSource, #unused#::Base.HasEltype, isz::Base.HasLength)
    @ Base ./array.jl:713
  [3] collect(itr::JUDI.MultiSource)
    @ Base ./array.jl:707
  [4] broadcastable(x::JUDI.MultiSource)
    @ Base.Broadcast ./broadcast.jl:717
  [5] broadcasted(::Function, ::JUDI.MultiSource)
    @ Base.Broadcast ./broadcast.jl:1309
  [6] envelope(x::judiVector{Float32, Matrix{Float32}}, y::judiVector{Float32, Matrix{Float32}})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:201
  [7] myloss(x::judiVector{Float32, Matrix{Float32}}, y::judiVector{Float32, Matrix{Float32}})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:205
  [8] objective_function(m_update::Matrix{Float64})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:251
  [9] (::SlimOptim.var"#objgrad!#11"{typeof(objective_function), result{Float32}})(g::Matrix{Float32}, x::Matrix{Float64})
    @ SlimOptim ~/Documents/Colada/r/julia-1.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:99
 [10] _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.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:172
 [11] 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.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:103
 [12] spg(funObj::Function, x::Matrix{Float32}, funProj::Function, options::SlimOptim.SPG_params)
    @ SlimOptim ~/Documents/Colada/r/julia-1.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:89
 [13] top-level scope
    @ ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:322
mloubout commented 3 months ago

Seems like need am few extra:

JUDI.make_input(jv::judiVector{T, Matrix{Complex{T}}}) where T = jv.data[1]

for the 1im .*

then just use norm(..., 2)^2 instead of sum(abs2) . If you really want sum(abs2 you will need to define abs2(x::judiMultiSourceVector (like abs)

kerim371 commented 3 months ago

Didn't work.

actually on the step of summing:

function JUDI.make_input(jv::judiVector{T, Matrix{Complex{T}}}) where T
    @info "JUDI.make_input"
    jv.data[1]
end

@everywhere function envelope(x, y)
    a = (x - y)
    b = 1im .* H(x - y)
    c = a .+ b  # hereI get error, RHS: judiVector{Float32, Matrix{ComplexF32}}, LHS: judiVector{Float32, Matrix{Float32}}
    return c
end

i get the error below.

And the function JUDI.make_input isn't called.

ERROR: InexactError: Float32(-3.3958447f-8 + 4.6234345f-9im)
Stacktrace:
  [1] Real
    @ ./complex.jl:44 [inlined]
  [2] convert
    @ ./number.jl:7 [inlined]
  [3] setindex!
    @ ./array.jl:969 [inlined]
  [4] _unsafe_copyto!(dest::Matrix{Float32}, doffs::Int64, src::Matrix{ComplexF32}, soffs::Int64, n::Int64)
    @ Base ./array.jl:250
  [5] unsafe_copyto!
    @ ./array.jl:304 [inlined]
  [6] _copyto_impl!
    @ ./array.jl:327 [inlined]
  [7] copyto!
    @ ./array.jl:314 [inlined]
  [8] copyto!
    @ ./array.jl:339 [inlined]
  [9] copyto!
    @ ./broadcast.jl:967 [inlined]
 [10] copyto!
    @ ./broadcast.jl:926 [inlined]
 [11] materialize!
    @ ./broadcast.jl:884 [inlined]
 [12] materialize!(dest::Matrix{Float32}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{2}, Nothing, typeof(identity), Tuple{Matrix{ComplexF32}}})
    @ Base.Broadcast ./broadcast.jl:881
 [13] materialize(bc::JUDI.MultiSource)
    @ JUDI ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/JEsVr/src/TimeModeling/Types/broadcasting.jl:46
 [14] envelope(x::judiVector{Float32, Matrix{Float32}}, y::judiVector{Float32, Matrix{Float32}})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:214
 [15] myloss(x::judiVector{Float32, Matrix{Float32}}, y::judiVector{Float32, Matrix{Float32}})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:223
 [16] objective_function(m_update::Matrix{Float64})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:269
 [17] (::SlimOptim.var"#objgrad!#11"{typeof(objective_function), result{Float32}})(g::Matrix{Float32}, x::Matrix{Float64})
    @ SlimOptim ~/Documents/Colada/r/julia-1.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:99
 [18] _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.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:172
 [19] 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.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:103
 [20] spg(funObj::Function, x::Matrix{Float32}, funProj::Function, options::SlimOptim.SPG_params)
    @ SlimOptim ~/Documents/Colada/r/julia-1.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:89
 [21] top-level scope
    @ ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:340
mloubout commented 3 months ago

What does c = b .+ a do ? Might be a small issue with type detection (i.e which is the higher type) I think it just picks the first one for the output type currently

kerim371 commented 3 months ago

What does c = b .+ a do ? Might be a small issue with type detection (i.e which is the higher type) I think it just picks the first one for the output type currently

It would be more correct to if I wrote something like that:

@everywhere function envelope(x, y)
    a = (x - y)
    b = 1im .* H(x - y)
    c = a .+ b  # here I get error, RHS: judiVector{Float32, Matrix{ComplexF32}}, LHS: judiVector{Float32, Matrix{Float32}}
    return norm(c, 2)^2
end

I splitted it to a and b to investigate the error (type of a and b and why JUDI.make_input in't called) .

Does JUDI.make_input should be called on the step b = 1im .* H(x - y) ? i.e. when multiplying by 1im.

For the now: H(x - y) gives the output of type judiVector{Float32, Matrix{ComplexF32}} 1im .* H(x - y) gives the output of type judiVector{Float32, Matrix{ComplexF32}} x - y gives the output of type judiVector{Float32, Matrix{Float32}}

mloubout commented 3 months ago

I mean.

 function envelope(x, y)
    a = (x - y)
    b = 1im .* H(x - y)
    c = b .+ a  # here I get error, RHS: judiVector{Float32, Matrix{ComplexF32}}, LHS: judiVector{Float32, Matrix{Float32}}
    return norm(c, 2)^2
end
kerim371 commented 3 months ago

Now seems that envelope function is computed fine but denvelop gives the error:

ERROR: InexactError: Float32(2.0317233f-5 - 5.650077f-10im)
Stacktrace:
  [1] Float32(z::ComplexF32)
    @ Base ./complex.jl:44
  [2] dot(a::judiVector{Float32, Matrix{ComplexF32}}, b::judiVector{Float32, Matrix{ComplexF32}})
    @ JUDI ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/JEsVr/src/TimeModeling/Types/abstract.jl:137
  [3] (::ChainRules.var"#1480#1484"{judiVector{Float32, Matrix{ComplexF32}}, judiVector{Float32, Matrix{ComplexF32}}, ChainRulesCore.ProjectTo{ComplexF64, NamedTuple{(), Tuple{}}}})()
    @ ChainRules ~/Documents/Colada/r/julia-1.9/.julia/packages/ChainRules/pEOSw/src/rulesets/Base/arraymath.jl:108
  [4] unthunk
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/ChainRulesCore/7MWx2/src/tangent_types/thunks.jl:204 [inlined]
  [5] wrap_chainrules_output
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/Zygote/WOy6z/src/compiler/chainrules.jl:110 [inlined]
  [6] map
    @ ./tuple.jl:275 [inlined]
  [7] wrap_chainrules_output
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/Zygote/WOy6z/src/compiler/chainrules.jl:111 [inlined]
  [8] ZBack
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/Zygote/WOy6z/src/compiler/chainrules.jl:211 [inlined]
  [9] (::Zygote.var"#3818#back#1211"{Zygote.ZBack{ChainRules.var"#times_pullback#1483"{Complex{Int64}, judiVector{Float32, Matrix{ComplexF32}}, ChainRulesCore.ProjectTo{AbstractArray, NamedTuple{(:element, :axes), Tuple{ChainRulesCore.ProjectTo{Float32, NamedTuple{(), Tuple{}}}, Tuple{Base.OneTo{Int64}}}}}, ChainRulesCore.ProjectTo{ComplexF64, NamedTuple{(), Tuple{}}}}}})(Δ::judiVector{Float32, Matrix{ComplexF32}})
    @ Zygote ~/Documents/Colada/r/julia-1.9/.julia/packages/ZygoteRules/M4xmc/src/adjoint.jl:72
 [10] Pullback
    @ ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:208 [inlined]
 [11] (::Zygote.Pullback{Tuple{typeof(envelope), judiVector{Float32, Matrix{Float32}}, judiVector{Float32, Matrix{Float32}}}, Tuple{Zygote.var"#3605#back#1090"{Zygote.var"#1086#1089"}, Zygote.Pullback{Tuple{typeof(H), judiVector{Float32, Matrix{Float32}}}, Tuple{Zygote.ZBack{Zygote.var"#plus_pullback#345"{Tuple{Float64, Int64}}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), judiVector{ComplexF32, Matrix{ComplexF32}}}, Tuple{}}, Zygote.ZBack{ChainRules.var"#/_pullback#1321"{Float64, Float64, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}}}, Zygote.ZBack{AbstractFFTsChainRulesCoreExt.var"#ifft_pullback#7"{Int64, ChainRulesCore.ProjectTo{AbstractArray, NamedTuple{(:element, :axes), Tuple{ChainRulesCore.ProjectTo{ComplexF32, NamedTuple{(), Tuple{}}}, Tuple{Base.OneTo{Int64}}}}}, Float32}}, Zygote.var"#4090#back#1347"{Zygote.var"#1343#1346"{Vector{Float64}}}, Zygote.var"#3802#back#1205"{Zygote.var"#1201#1204"{Vector{Float32}, judiVector{ComplexF32, Matrix{ComplexF32}}}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.broadcasted), typeof(sign), StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Tuple{Zygote.var"#2173#back#293"{Zygote.var"#291#292"{Tuple{Tuple{Nothing, Nothing, Nothing}, Tuple{}}, Zygote.var"#4123#back#1360"{Zygote.var"#bc_fwd_back#1398"{Vector{ForwardDiff.Dual{Nothing, Float64, 1}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Val{1}}}}}, Zygote.var"#2017#back#204"{typeof(identity)}, Zygote.var"#2017#back#204"{typeof(identity)}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.broadcastable), StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Tuple{}}, Zygote.var"#2173#back#293"{Zygote.var"#291#292"{Tuple{Tuple{Nothing}, Tuple{}}, Zygote.var"#combine_styles_pullback#1166"{Tuple{Nothing, Nothing}}}}, Zygote.var"#2849#back#672"{Zygote.var"#map_back#666"{typeof(Base.Broadcast.broadcastable), 1, Tuple{Tuple{}}, Tuple{Val{0}}, Tuple{}}}}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), Vector{Float32}}, Tuple{}}, Zygote.ZBack{ChainRules.var"#size_pullback#919"}, Zygote.ZBack{AbstractFFTsChainRulesCoreExt.var"#fft_pullback#1"{Int64, ChainRulesCore.ProjectTo{AbstractArray, NamedTuple{(:element, :axes), Tuple{ChainRulesCore.ProjectTo{Float32, NamedTuple{(), Tuple{}}}, Tuple{Base.OneTo{Int64}}}}}}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), Vector{Float64}}, Tuple{}}, Zygote.var"#3059#back#801"{Zygote.var"#797#800"}, Zygote.ZBack{ChainRules.var"#-_pullback#1329"{Int64, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}}}, Zygote.Pullback{Tuple{typeof(ifftshift), Vector{Float64}}, Tuple{Zygote.Pullback{Tuple{typeof(ndims), Vector{Float64}}, Tuple{}}, Zygote.ZBack{AbstractFFTsChainRulesCoreExt.var"#ifftshift_pullback#20"{UnitRange{Int64}, ChainRulesCore.ProjectTo{AbstractArray, NamedTuple{(:element, :axes), Tuple{ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}, Tuple{Base.OneTo{Int64}}}}}}}, Zygote.ZBack{ChainRules.var"#:_pullback#276"{Tuple{Int64, Int64}}}}}, Zygote.ZBack{ChainRules.var"#/_pullback#1321"{Float64, Float64, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}}}, Zygote.ZBack{ChainRules.var"#:_pullback#276"{Tuple{Float64, Float64}}}}}, Zygote.var"#3605#back#1090"{Zygote.var"#1086#1089"}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), judiVector{Float32, Matrix{ComplexF32}}}, Tuple{}}, Zygote.var"#3754#back#1181"{Zygote.var"#1175#1179"{Tuple{judiVector{Float32, Matrix{ComplexF32}}, judiVector{Float32, Matrix{Float32}}}}}, Zygote.var"#1990#back#194"{Zygote.var"#190#193"{Zygote.Context{false}, GlobalRef, Complex{Bool}}}, Zygote.ZBack{ChainRules.var"#norm_pullback_p#1999"{judiVector{Float32, Matrix{ComplexF32}}, Int64, Float32}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), judiVector{Float32, Matrix{ComplexF32}}}, Tuple{}}, Zygote.var"#1926#back#161"{Zygote.var"#157#160"}, Zygote.ZBack{Zygote.var"#literal_pow_pullback#331"{2, Float32}}, Zygote.ZBack{ChainRules.var"#times_pullback2#1331"{Int64, Complex{Bool}}}, Zygote.var"#3818#back#1211"{Zygote.ZBack{ChainRules.var"#times_pullback#1483"{Complex{Int64}, judiVector{Float32, Matrix{ComplexF32}}, ChainRulesCore.ProjectTo{AbstractArray, NamedTuple{(:element, :axes), Tuple{ChainRulesCore.ProjectTo{Float32, NamedTuple{(), Tuple{}}}, Tuple{Base.OneTo{Int64}}}}}, ChainRulesCore.ProjectTo{ComplexF64, NamedTuple{(), Tuple{}}}}}}}})(Δ::Float32)
    @ Zygote ~/Documents/Colada/r/julia-1.9/.julia/packages/Zygote/WOy6z/src/compiler/interface2.jl:0
 [12] Pullback
    @ ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:213 [inlined]
 [13] (::Zygote.Pullback{Tuple{var"#17#18"{judiVector{Float32, Matrix{Float32}}}, judiVector{Float32, Matrix{Float32}}}, Tuple{Zygote.var"#2184#back#303"{Zygote.var"#back#302"{:y, Zygote.Context{false}, var"#17#18"{judiVector{Float32, Matrix{Float32}}}, judiVector{Float32, Matrix{Float32}}}}, Zygote.Pullback{Tuple{typeof(envelope), judiVector{Float32, Matrix{Float32}}, judiVector{Float32, Matrix{Float32}}}, Tuple{Zygote.var"#3605#back#1090"{Zygote.var"#1086#1089"}, Zygote.Pullback{Tuple{typeof(H), judiVector{Float32, Matrix{Float32}}}, Tuple{Zygote.ZBack{Zygote.var"#plus_pullback#345"{Tuple{Float64, Int64}}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), judiVector{ComplexF32, Matrix{ComplexF32}}}, Tuple{}}, Zygote.ZBack{ChainRules.var"#/_pullback#1321"{Float64, Float64, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}}}, Zygote.ZBack{AbstractFFTsChainRulesCoreExt.var"#ifft_pullback#7"{Int64, ChainRulesCore.ProjectTo{AbstractArray, NamedTuple{(:element, :axes), Tuple{ChainRulesCore.ProjectTo{ComplexF32, NamedTuple{(), Tuple{}}}, Tuple{Base.OneTo{Int64}}}}}, Float32}}, Zygote.var"#4090#back#1347"{Zygote.var"#1343#1346"{Vector{Float64}}}, Zygote.var"#3802#back#1205"{Zygote.var"#1201#1204"{Vector{Float32}, judiVector{ComplexF32, Matrix{ComplexF32}}}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.broadcasted), typeof(sign), StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Tuple{Zygote.var"#2173#back#293"{Zygote.var"#291#292"{Tuple{Tuple{Nothing, Nothing, Nothing}, Tuple{}}, Zygote.var"#4123#back#1360"{Zygote.var"#bc_fwd_back#1398"{Vector{ForwardDiff.Dual{Nothing, Float64, 1}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Val{1}}}}}, Zygote.var"#2017#back#204"{typeof(identity)}, Zygote.var"#2017#back#204"{typeof(identity)}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.broadcastable), StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Tuple{}}, Zygote.var"#2173#back#293"{Zygote.var"#291#292"{Tuple{Tuple{Nothing}, Tuple{}}, Zygote.var"#combine_styles_pullback#1166"{Tuple{Nothing, Nothing}}}}, Zygote.var"#2849#back#672"{Zygote.var"#map_back#666"{typeof(Base.Broadcast.broadcastable), 1, Tuple{Tuple{}}, Tuple{Val{0}}, Tuple{}}}}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), Vector{Float32}}, Tuple{}}, Zygote.ZBack{ChainRules.var"#size_pullback#919"}, Zygote.ZBack{AbstractFFTsChainRulesCoreExt.var"#fft_pullback#1"{Int64, ChainRulesCore.ProjectTo{AbstractArray, NamedTuple{(:element, :axes), Tuple{ChainRulesCore.ProjectTo{Float32, NamedTuple{(), Tuple{}}}, Tuple{Base.OneTo{Int64}}}}}}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), Vector{Float64}}, Tuple{}}, Zygote.var"#3059#back#801"{Zygote.var"#797#800"}, Zygote.ZBack{ChainRules.var"#-_pullback#1329"{Int64, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}}}, Zygote.Pullback{Tuple{typeof(ifftshift), Vector{Float64}}, Tuple{Zygote.Pullback{Tuple{typeof(ndims), Vector{Float64}}, Tuple{}}, Zygote.ZBack{AbstractFFTsChainRulesCoreExt.var"#ifftshift_pullback#20"{UnitRange{Int64}, ChainRulesCore.ProjectTo{AbstractArray, NamedTuple{(:element, :axes), Tuple{ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}, Tuple{Base.OneTo{Int64}}}}}}}, Zygote.ZBack{ChainRules.var"#:_pullback#276"{Tuple{Int64, Int64}}}}}, Zygote.ZBack{ChainRules.var"#/_pullback#1321"{Float64, Float64, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}}}, Zygote.ZBack{ChainRules.var"#:_pullback#276"{Tuple{Float64, Float64}}}}}, Zygote.var"#3605#back#1090"{Zygote.var"#1086#1089"}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), judiVector{Float32, Matrix{ComplexF32}}}, Tuple{}}, Zygote.var"#3754#back#1181"{Zygote.var"#1175#1179"{Tuple{judiVector{Float32, Matrix{ComplexF32}}, judiVector{Float32, Matrix{Float32}}}}}, Zygote.var"#1990#back#194"{Zygote.var"#190#193"{Zygote.Context{false}, GlobalRef, Complex{Bool}}}, Zygote.ZBack{ChainRules.var"#norm_pullback_p#1999"{judiVector{Float32, Matrix{ComplexF32}}, Int64, Float32}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), judiVector{Float32, Matrix{ComplexF32}}}, Tuple{}}, Zygote.var"#1926#back#161"{Zygote.var"#157#160"}, Zygote.ZBack{Zygote.var"#literal_pow_pullback#331"{2, Float32}}, Zygote.ZBack{ChainRules.var"#times_pullback2#1331"{Int64, Complex{Bool}}}, Zygote.var"#3818#back#1211"{Zygote.ZBack{ChainRules.var"#times_pullback#1483"{Complex{Int64}, judiVector{Float32, Matrix{ComplexF32}}, ChainRulesCore.ProjectTo{AbstractArray, NamedTuple{(:element, :axes), Tuple{ChainRulesCore.ProjectTo{Float32, NamedTuple{(), Tuple{}}}, Tuple{Base.OneTo{Int64}}}}}, ChainRulesCore.ProjectTo{ComplexF64, NamedTuple{(), Tuple{}}}}}}}}}})(Δ::Float32)
    @ Zygote ~/Documents/Colada/r/julia-1.9/.julia/packages/Zygote/WOy6z/src/compiler/interface2.jl:0
 [14] (::Zygote.var"#75#76"{Zygote.Pullback{Tuple{var"#17#18"{judiVector{Float32, Matrix{Float32}}}, judiVector{Float32, Matrix{Float32}}}, Tuple{Zygote.var"#2184#back#303"{Zygote.var"#back#302"{:y, Zygote.Context{false}, var"#17#18"{judiVector{Float32, Matrix{Float32}}}, judiVector{Float32, Matrix{Float32}}}}, Zygote.Pullback{Tuple{typeof(envelope), judiVector{Float32, Matrix{Float32}}, judiVector{Float32, Matrix{Float32}}}, Tuple{Zygote.var"#3605#back#1090"{Zygote.var"#1086#1089"}, Zygote.Pullback{Tuple{typeof(H), judiVector{Float32, Matrix{Float32}}}, Tuple{Zygote.ZBack{Zygote.var"#plus_pullback#345"{Tuple{Float64, Int64}}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), judiVector{ComplexF32, Matrix{ComplexF32}}}, Tuple{}}, Zygote.ZBack{ChainRules.var"#/_pullback#1321"{Float64, Float64, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}}}, Zygote.ZBack{AbstractFFTsChainRulesCoreExt.var"#ifft_pullback#7"{Int64, ChainRulesCore.ProjectTo{AbstractArray, NamedTuple{(:element, :axes), Tuple{ChainRulesCore.ProjectTo{ComplexF32, NamedTuple{(), Tuple{}}}, Tuple{Base.OneTo{Int64}}}}}, Float32}}, Zygote.var"#4090#back#1347"{Zygote.var"#1343#1346"{Vector{Float64}}}, Zygote.var"#3802#back#1205"{Zygote.var"#1201#1204"{Vector{Float32}, judiVector{ComplexF32, Matrix{ComplexF32}}}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.broadcasted), typeof(sign), StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Tuple{Zygote.var"#2173#back#293"{Zygote.var"#291#292"{Tuple{Tuple{Nothing, Nothing, Nothing}, Tuple{}}, Zygote.var"#4123#back#1360"{Zygote.var"#bc_fwd_back#1398"{Vector{ForwardDiff.Dual{Nothing, Float64, 1}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Val{1}}}}}, Zygote.var"#2017#back#204"{typeof(identity)}, Zygote.var"#2017#back#204"{typeof(identity)}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.broadcastable), StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Tuple{}}, Zygote.var"#2173#back#293"{Zygote.var"#291#292"{Tuple{Tuple{Nothing}, Tuple{}}, Zygote.var"#combine_styles_pullback#1166"{Tuple{Nothing, Nothing}}}}, Zygote.var"#2849#back#672"{Zygote.var"#map_back#666"{typeof(Base.Broadcast.broadcastable), 1, Tuple{Tuple{}}, Tuple{Val{0}}, Tuple{}}}}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), Vector{Float32}}, Tuple{}}, Zygote.ZBack{ChainRules.var"#size_pullback#919"}, Zygote.ZBack{AbstractFFTsChainRulesCoreExt.var"#fft_pullback#1"{Int64, ChainRulesCore.ProjectTo{AbstractArray, NamedTuple{(:element, :axes), Tuple{ChainRulesCore.ProjectTo{Float32, NamedTuple{(), Tuple{}}}, Tuple{Base.OneTo{Int64}}}}}}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), Vector{Float64}}, Tuple{}}, Zygote.var"#3059#back#801"{Zygote.var"#797#800"}, Zygote.ZBack{ChainRules.var"#-_pullback#1329"{Int64, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}}}, Zygote.Pullback{Tuple{typeof(ifftshift), Vector{Float64}}, Tuple{Zygote.Pullback{Tuple{typeof(ndims), Vector{Float64}}, Tuple{}}, Zygote.ZBack{AbstractFFTsChainRulesCoreExt.var"#ifftshift_pullback#20"{UnitRange{Int64}, ChainRulesCore.ProjectTo{AbstractArray, NamedTuple{(:element, :axes), Tuple{ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}, Tuple{Base.OneTo{Int64}}}}}}}, Zygote.ZBack{ChainRules.var"#:_pullback#276"{Tuple{Int64, Int64}}}}}, Zygote.ZBack{ChainRules.var"#/_pullback#1321"{Float64, Float64, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}, ChainRulesCore.ProjectTo{Float64, NamedTuple{(), Tuple{}}}}}, Zygote.ZBack{ChainRules.var"#:_pullback#276"{Tuple{Float64, Float64}}}}}, Zygote.var"#3605#back#1090"{Zygote.var"#1086#1089"}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), judiVector{Float32, Matrix{ComplexF32}}}, Tuple{}}, Zygote.var"#3754#back#1181"{Zygote.var"#1175#1179"{Tuple{judiVector{Float32, Matrix{ComplexF32}}, judiVector{Float32, Matrix{Float32}}}}}, Zygote.var"#1990#back#194"{Zygote.var"#190#193"{Zygote.Context{false}, GlobalRef, Complex{Bool}}}, Zygote.ZBack{ChainRules.var"#norm_pullback_p#1999"{judiVector{Float32, Matrix{ComplexF32}}, Int64, Float32}}, Zygote.Pullback{Tuple{typeof(Base.Broadcast.materialize), judiVector{Float32, Matrix{ComplexF32}}}, Tuple{}}, Zygote.var"#1926#back#161"{Zygote.var"#157#160"}, Zygote.ZBack{Zygote.var"#literal_pow_pullback#331"{2, Float32}}, Zygote.ZBack{ChainRules.var"#times_pullback2#1331"{Int64, Complex{Bool}}}, Zygote.var"#3818#back#1211"{Zygote.ZBack{ChainRules.var"#times_pullback#1483"{Complex{Int64}, judiVector{Float32, Matrix{ComplexF32}}, ChainRulesCore.ProjectTo{AbstractArray, NamedTuple{(:element, :axes), Tuple{ChainRulesCore.ProjectTo{Float32, NamedTuple{(), Tuple{}}}, Tuple{Base.OneTo{Int64}}}}}, ChainRulesCore.ProjectTo{ComplexF64, NamedTuple{(), Tuple{}}}}}}}}}}})(Δ::Float32)
    @ Zygote ~/Documents/Colada/r/julia-1.9/.julia/packages/Zygote/WOy6z/src/compiler/interface.jl:45
 [15] gradient(f::Function, args::judiVector{Float32, Matrix{Float32}})
    @ Zygote ~/Documents/Colada/r/julia-1.9/.julia/packages/Zygote/WOy6z/src/compiler/interface.jl:97
 [16] denvelope(x::judiVector{Float32, Matrix{Float32}}, y::judiVector{Float32, Matrix{Float32}})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:213
 [17] myloss(x::judiVector{Float32, Matrix{Float32}}, y::judiVector{Float32, Matrix{Float32}})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:215
 [18] objective_function(m_update::Matrix{Float64})
    @ Main ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:261
 [19] (::SlimOptim.var"#objgrad!#11"{typeof(objective_function), result{Float32}})(g::Matrix{Float32}, x::Matrix{Float64})
    @ SlimOptim ~/Documents/Colada/r/julia-1.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:99
 [20] _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.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:172
 [21] 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.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:103
 [22] spg(funObj::Function, x::Matrix{Float32}, funProj::Function, options::SlimOptim.SPG_params)
    @ SlimOptim ~/Documents/Colada/r/julia-1.9/.julia/packages/SlimOptim/33Fm9/src/SPGSlim.jl:89
 [23] top-level scope
    @ ~/Documents/JUDI.jl/examples/field_examples/viking_graben_line12/fwi/fwi_spg.jl:332
mloubout commented 3 months ago

https://github.com/slimgroup/JUDI.jl/blob/05e9ff96800426db0988c231e2366bff9075414a/src/TimeModeling/Types/abstract.jl#L137

That's your issue since T is Float32, probably need T(real) (assuming dot returns a real)

kerim371 commented 3 months ago

ok, I will try to investigate that