search

slimgroup / JUDI.jl

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

New version `subsampling_factor` doesn't work #242

Closed kerim371 closed 5 months ago

kerim371 commented 5 months ago

Hi,

I'm testing new version of JUDI and it seems that subsampling_factor=16 gives an error during FWI:

ERROR: PyError ($(Expr(:escape, :(ccall(#= /home/kerim/Documents/Colada/r/julia-1.9/.julia/packages/PyCall/1gn3u/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))) <class 'ValueError'>
ValueError('Override `us_u(t_sub, x, y)` is incompatible with overrides `[damp(x, y), m(x, y), rcvu(time, p_rcvu), srcu(time, p_srcu), u(t, x, y), us_u(t_sub, x, y)]`')
  File "/home/kerim/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/pysource/interface.py", line 359, in J_adjoint
    return J_adjoint_standard(model, src_coords, wavelet, rec_coords, recin,
  File "/home/kerim/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/pysource/interface.py", line 478, in J_adjoint_standard
    rec, u, Iu, _ = ffunc(model, src_coords, rec_coords, wavelet, save=True, nlind=nlind,
  File "/home/kerim/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/pysource/propagators.py", line 77, in forward
    summary = op(**kw)
  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/operator/operator.py", line 753, in __call__
    return self.apply(**kwargs)
  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/operator/operator.py", line 819, in apply
    args = self.arguments(**kwargs)
  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/operator/operator.py", line 660, in arguments
    args = self._prepare_arguments(**kwargs)
  File "/home/kerim/Documents/Colada/r/python-install/lib/python3.9/site-packages/devito/operator/operator.py", line 543, in _prepare_arguments
    raise ValueError("Override `%s` is incompatible with overrides `%s`" %

Stacktrace:
  [1] pyerr_check
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/PyCall/1gn3u/src/exception.jl:75 [inlined]
  [2] pyerr_check
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/PyCall/1gn3u/src/exception.jl:79 [inlined]
  [3] _handle_error(msg::String)
    @ PyCall ~/Documents/Colada/r/julia-1.9/.julia/packages/PyCall/1gn3u/src/exception.jl:96
  [4] macro expansion
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/PyCall/1gn3u/src/exception.jl:110 [inlined]
  [5] #107
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/PyCall/1gn3u/src/pyfncall.jl:43 [inlined]
  [6] disable_sigint
    @ ./c.jl:473 [inlined]
  [7] __pycall!
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/PyCall/1gn3u/src/pyfncall.jl:42 [inlined]
  [8] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Vararg{Matrix{Float32}, 4}}, nargs::Int64, kw::PyCall.PyObject)
    @ PyCall ~/Documents/Colada/r/julia-1.9/.julia/packages/PyCall/1gn3u/src/pyfncall.jl:29
  [9] _pycall!(ret::PyCall.PyObject, o::PyCall.PyObject, args::Tuple{PyCall.PyObject, Vararg{Matrix{Float32}, 4}}, kwargs::Base.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.9/.julia/packages/PyCall/1gn3u/src/pyfncall.jl:11
 [10] #pycall#112
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/PyCall/1gn3u/src/pyfncall.jl:80 [inlined]
 [11] pycall
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/PyCall/1gn3u/src/pyfncall.jl:80 [inlined]
 [12] #4
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/JUDI.jl:82 [inlined]
 [13] (::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.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, Vararg{Matrix{Float32}, 4}}}})()
    @ JUDI ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/JUDI.jl:74
 [14] lock(f::JUDI.var"#1#2"{JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.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, Vararg{Matrix{Float32}, 4}}}}, l::ReentrantLock)
    @ Base ./lock.jl:229
 [15] pylock(f::JUDI.var"#4#5"{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}, Base.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, Vararg{Matrix{Float32}, 4}}})
    @ JUDI ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/JUDI.jl:71
 [16] rlock_pycall(::PyCall.PyObject, ::Type{Tuple{Float32, PyCall.PyArray, PyCall.PyObject, PyCall.PyObject}}, ::PyCall.PyObject, ::Vararg{Any}; kw::Base.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.9/.julia/packages/JUDI/9qfzJ/src/JUDI.jl:81
 [17] macro expansion
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/TimeModeling/Modeling/misfit_fg.jl:50 [inlined]
 [18] macro expansion
    @ ./timing.jl:393 [inlined]
 [19] macro expansion
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/JUDI.jl:141 [inlined]
 [20] 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.9/.julia/packages/JUDI/9qfzJ/src/TimeModeling/Modeling/misfit_fg.jl:49
 [21] 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)
    @ JUDI ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/TimeModeling/Modeling/misfit_fg.jl:78
 [22] macro expansion
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/TimeModeling/Modeling/propagation.jl:36 [inlined]
 [23] macro expansion
    @ ./timing.jl:393 [inlined]
 [24] macro expansion
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/JUDI.jl:141 [inlined]
 [25] run_and_reduce(func::Function, #unused#::Nothing, nsrc::Int64, arg_func::JUDI.var"#276#277"{JUDIOptions, Base.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:nlind, :lin), Tuple{Bool, Bool}}}, JUDI.IsoModel{Float32, 2}, judiVector{Float32, Matrix{Float32}}, judiVector{Float32, Matrix{Float32}}, Nothing})
    @ JUDI ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/TimeModeling/Modeling/propagation.jl:35
 [26] 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.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:nlind, :lin), Tuple{Bool, Bool}}})
    @ JUDI ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/TimeModeling/Modeling/propagation.jl:108
 [27] multi_src_fg!
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/TimeModeling/Modeling/propagation.jl:100 [inlined]
 [28] #multi_exp_fg!#255
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/TimeModeling/Modeling/misfit_fg.jl:191 [inlined]
 [29] multi_exp_fg!
    @ ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/TimeModeling/Modeling/misfit_fg.jl:191 [inlined]
 [30] 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.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ JUDI ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/TimeModeling/Modeling/misfit_fg.jl:172
 [31] fwi_objective(model::JUDI.IsoModel{Float32, 2}, q::judiVector{Float32, Matrix{Float32}}, dobs::judiVector{Float32, Matrix{Float32}}; options::JUDIOptions, kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ JUDI ~/Documents/Colada/r/julia-1.9/.julia/packages/JUDI/9qfzJ/src/TimeModeling/Modeling/misfit_fg.jl:137
 [32] top-level scope
    @ ~/Documents/Colada/r/julia-1.9/.julia/dev/JUDI/examples/scripts/fwi_example_NLopt.jl:63

The script that I run is modified fwi_example_NLopt.jl (function f!(x,grad) is intensionally commited so that julia outputs the error message):

# 2D FWI on Overthrust model with L-BFGS using NLopt library
# Author: Philipp Witte, pwitte@eoas.ubc.ca
# Date: December 2017
#

using Statistics, Random, LinearAlgebra
using JUDI, HDF5, NLopt, SegyIO

# Load starting model
n,d,o,m0 = read(h5open("$(JUDI.JUDI_DATA)/overthrust_model.h5","r"), "n", "d", "o", "m0")
model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)

# Bound constraints
v0 = sqrt.(1f0 ./m0)
vmin = ones(Float32, model0.n) .* 1.3f0
vmax = ones(Float32, model0.n) .* 6.5f0

# Slowness squared [s^2/km^2]
mmin = vec((1f0 ./ vmax).^2)
mmax = vec((1f0 ./ vmin).^2)

# Load data
block = segy_read("$(JUDI.JUDI_DATA)/overthrust_shot_records.segy")
d_obs = judiVector(block)

# Set up wavelet
src_geometry = Geometry(block; key="source")
wavelet = ricker_wavelet(src_geometry.t[1], src_geometry.dt[1], 0.008f0)    # 8 Hz wavelet
q = judiVector(src_geometry, wavelet)

# t0 somehow can't be 0
Ml_ref = judiDataMute(q.geometry, d_obs.geometry, vp=1100f0, t0=0.001f0, mode=:reflection) # keep reflections
Ml_tur = judiDataMute(q.geometry, d_obs.geometry, vp=1100f0, t0=0.001f0, mode=:turning)    # keep turning waves

# JUDI options
global jopt = JUDI.Options(
    IC = "fwi",
    # limit_m = true,
    # buffer_size = 1000,
    # optimal_checkpointing=false,
    subsampling_factor=16,
    # free_surface=true,
    # space_order=32
    )

############################### FWI ###########################################

# optimization parameters
batchsize = 16
count = 0

# make sure that commited function `f!` will print error as well as the script fails
x = copy(model0.m)
grad = zeros(Float32, model0.n)

# NLopt objective function
println("No.  ", "fval         ", "norm(gradient)")
# function f!(x,grad)

    # Update model
    model0.m .= convert(Array{Float32, 2}, reshape(x, model0.n))

    # Seclect batch and calculate gradient
    i = randperm(d_obs.nsrc)[1:batchsize]
    fval, gradient = fwi_objective(model0, q[i], Ml_ref[i]*d_obs[i], options=jopt)

    # Reset gradient in water column to zero
    gradient = reshape(gradient, model0.n)
    gradient[:, 1:21] .= 0f0
    grad[1:end] = vec(gradient)

    global count; count += 1
    println(count, "    ", fval, "    ", norm(grad))
    return convert(Float64, fval)
# end

# Optimization parameters
opt = Opt(:LD_LBFGS, prod(model0.n))
lower_bounds!(opt, mmin); upper_bounds!(opt, mmax)
min_objective!(opt, f!)
maxeval!(opt, parse(Int, get(ENV, "NITER", "10")))
(minf, minx, ret) = optimize(opt, vec(model0.m.data))

JUDI: v3.3.10

mloubout commented 5 months ago

I can't reproduce it it runs fine for me. What version of devito are you using?

kerim371 commented 5 months ago

@mloubout right!

I've been using devito v4.8.2 and it used to give error. But after I updated to devito v4.8.3 it started to work fine.

Thank you!