From 350622c81ea07780e7fec2cef377614f27de904a Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 23 Oct 2024 12:21:03 -0400 Subject: [PATCH] fix dm input --- .github/workflows/ci-examples.yml | 9 ++-- .github/workflows/ci-judi.yml | 9 ++-- .github/workflows/ci-op.yml | 9 ++-- examples/scripts/modeling_basic_2D.jl | 16 +++++-- src/TimeModeling/Modeling/losses.jl | 22 ++------- src/TimeModeling/Modeling/misfit_fg.jl | 3 -- src/TimeModeling/Modeling/python_interface.jl | 24 ++++------ src/TimeModeling/Modeling/twri_objective.jl | 7 --- .../Preconditioners/DataPreconditioners.jl | 10 ++--- src/TimeModeling/Preconditioners/base.jl | 1 + src/TimeModeling/Utils/auxiliaryFunctions.jl | 18 +++----- src/TimeModeling/Utils/time_utilities.jl | 2 + src/pysource/fields.py | 9 ++-- src/pysource/geom_utils.py | 3 +- src/pysource/interface.py | 7 ++- src/pysource/models.py | 27 +++++++---- src/pysource/propagators.py | 2 +- src/pysource/sensitivity.py | 8 ++-- src/pysource/utils.py | 8 ++++ src/utilities.jl | 4 +- test/test_basics.jl | 45 ++++++++++++------- test/test_preconditioners.jl | 6 ++- 22 files changed, 136 insertions(+), 113 deletions(-) diff --git a/.github/workflows/ci-examples.yml b/.github/workflows/ci-examples.yml index 35a294e90..cdc82c75c 100644 --- a/.github/workflows/ci-examples.yml +++ b/.github/workflows/ci-examples.yml @@ -42,6 +42,7 @@ jobs: DEVITO_LANGUAGE: "openmp" DEVITO_LOGGING: "ERROR" OMP_NUM_THREADS: 1 + JULIA_NUM_THREADS: 1 NITER: 2 strategy: @@ -60,21 +61,23 @@ jobs: version: ${{ matrix.version }} arch: x64 + - uses: julia-actions/cache@v2 + - name: Set julia python run: | python3 -m pip install -U pip python3 -m pip install "matplotlib<3.9" seiscm colorcet echo "JULIA_PYTHONCALL_EXE=$(which python3)" >> $GITHUB_ENV echo "JULIA_CONDAPKG_BACKEND=\"Null\"" >> $GITHUB_ENV - python3 -m pip install devito[tests,extras]@git+https://github.com/devitocodes/devito.git - JULIA_PYTHONCALL_EXE=$(which python3) JULIA_CONDAPKG_BACKEND="Null" julia -e 'using Pkg;Pkg.add(["PythonCall", "PyPlot", "SlimPlotting"]);Pkg.build("PythonCall")' + echo "PYTHON=$(which python3)" >> $GITHUB_ENV + echo "PYCALL_JL_RUNTIME_PYTHON=$(which python3)" >> $GITHUB_ENV - name: Build JUDI uses: julia-actions/julia-buildpkg@latest - name: Install packages run: | - julia -e 'using Pkg;Pkg.add(["NLopt", "Flux", "JOLI", "Zygote", "IterativeSolvers", "SlimOptim", "HDF5", "SegyIO", "SetIntersectionProjection"])' + julia -e 'using Pkg;Pkg.add(["SlimPlotting", "PyPlot", "NLopt", "Flux", "JOLI", "Zygote", "IterativeSolvers", "SlimOptim", "HDF5", "SegyIO", "SetIntersectionProjection"])' julia -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd()))' - name: Run examples diff --git a/.github/workflows/ci-judi.yml b/.github/workflows/ci-judi.yml index 4a98c46c9..9b41553b2 100644 --- a/.github/workflows/ci-judi.yml +++ b/.github/workflows/ci-judi.yml @@ -22,6 +22,7 @@ jobs: DEVITO_ARCH: ${{ matrix.cc }} DEVITO_LANGUAGE: "openmp" OMP_NUM_THREADS: 4 + JULIA_NUM_THREADS: 1 GROUP: "JUDI" strategy: @@ -49,11 +50,13 @@ jobs: version: ${{ matrix.version }} arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v2 + - name: Setup clang for osx if: runner.os == 'macOS' run: | brew install llvm libomp - echo "/opt/homebrew/bin:/opt/homebrew/opt/llvm/bin" >> $GITHUB_PATH + echo "/opt/homebrew/opt/llvm/bin" >> $GITHUB_PATH - name: Set up Python 3.9 uses: actions/setup-python@v5 @@ -62,10 +65,10 @@ jobs: - name: Set julia python run: | + echo "PYTHON=$(which python3)" >> $GITHUB_ENV + echo "PYCALL_JL_RUNTIME_PYTHON=$(which python3)" >> $GITHUB_ENV echo "JULIA_PYTHONCALL_EXE=$(which python3)" >> $GITHUB_ENV echo "JULIA_CONDAPKG_BACKEND=\"Null\"" >> $GITHUB_ENV - python3 -m pip install devito[tests,extras]@git+https://github.com/devitocodes/devito.git - JULIA_PYTHONCALL_EXE=$(which python3) JULIA_CONDAPKG_BACKEND="Null" julia -e 'using Pkg;Pkg.add("PythonCall");Pkg.build("PythonCall")' - name: Build JUDI uses: julia-actions/julia-buildpkg@latest diff --git a/.github/workflows/ci-op.yml b/.github/workflows/ci-op.yml index 875f72122..db4a82c33 100644 --- a/.github/workflows/ci-op.yml +++ b/.github/workflows/ci-op.yml @@ -23,6 +23,7 @@ jobs: DEVITO_LANGUAGE: "openmp" DEVITO_LOGGING: "INFO" OMP_NUM_THREADS: ${{ matrix.omp }} + JULIA_NUM_THREADS: 1 GROUP: ${{ matrix.op }} strategy: @@ -82,11 +83,13 @@ jobs: version: ${{ matrix.version }} arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v2 + - name: Setup clang for osx if: runner.os == 'macOS' run: | brew install llvm libomp - echo "/opt/homebrew/bin:/opt/homebrew/opt/llvm/bin" >> $GITHUB_PATH + echo "/opt/homebrew/opt/llvm/bin" >> $GITHUB_PATH - name: Set up Python 3.9 uses: actions/setup-python@v5 @@ -95,10 +98,10 @@ jobs: - name: Set julia python run: | + echo "PYTHON=$(which python3)" >> $GITHUB_ENV + echo "PYCALL_JL_RUNTIME_PYTHON=$(which python3)" >> $GITHUB_ENV echo "JULIA_PYTHONCALL_EXE=$(which python3)" >> $GITHUB_ENV echo "JULIA_CONDAPKG_BACKEND=\"Null\"" >> $GITHUB_ENV - python3 -m pip install devito[tests,extras]@git+https://github.com/devitocodes/devito.git - JULIA_PYTHONCALL_EXE=$(which python3) JULIA_CONDAPKG_BACKEND="Null" julia -e 'using Pkg;Pkg.add("PythonCall");Pkg.build("PythonCall")' - name: Build JUDI uses: julia-actions/julia-buildpkg@latest diff --git a/examples/scripts/modeling_basic_2D.jl b/examples/scripts/modeling_basic_2D.jl index 99d98bf1d..c6beb6f6f 100644 --- a/examples/scripts/modeling_basic_2D.jl +++ b/examples/scripts/modeling_basic_2D.jl @@ -137,7 +137,7 @@ dot1 = dot(q, qad) # dot2 = dot(dobs, dobs) # Compare -@show dot1, dot2, (dot2 - dot2)/(dot1 + dot2) +@show dot1, dot2, (dot1 - dot2)/(dot1 + dot2) #' # Inversion #' Our main goal is to provide an inversion framework for seismic inversion. To this end, as shown earlier, @@ -151,15 +151,25 @@ rtm = adjoint(J)*dD #' We show the linearized data. fig = figure() -plot_sdata(dobs[2]; new_fig=false, name="Linearized data", cmap=dcmap) +plot_sdata(dD[2]; new_fig=false, name="Linearized data", cmap=dcmap) display(fig) - #' And the RTM image fig = figure() plot_simage(rtm'; new_fig=false, name="RTM image", cmap=imcmap) display(fig) +#' We can easily now again test the adjointness of our operator with the standard dot test. Because we +#' intend to conserve our linear algebra abstraction, `judiVector` implements all the necessary linear +#' algebra functions such as dot product or norm to be used directly. +# +dot3 = dot(dm, rtm) +# +dot4 = dot(dD, dD) +# Compare +@show dot3, dot4, (dot3 - dot4)/(dot3 + dot4) + + #' ## Inversion utility functions #' We currently introduced the lineaar operators that allow to write seismic modeling and inversion in a high-level, linear algebra way. These linear operators allow the script to closely follow the mathematics and to be readable and understandable. #' diff --git a/src/TimeModeling/Modeling/losses.jl b/src/TimeModeling/Modeling/losses.jl index 376543306..d95a40d24 100644 --- a/src/TimeModeling/Modeling/losses.jl +++ b/src/TimeModeling/Modeling/losses.jl @@ -12,19 +12,12 @@ and its derivative w.r.t `x` `x-y` """ -function _mse(x::AT, y::AbstractArray{T}) where {T<:Number, AT<:AbstractArray{T}} - f = .5f0 * norm(x - y, 2)^2 +function mse(x::AbstractArray{T}, y::AbstractArray{T}) where {T<:Number} + f = T(.5) * norm(x - y, 2)^2 r = x - y return f, r end -function mse(x::PyArray{T}, y::AbstractArray{T}) where {T<:Number} - f, r = _mse(x, y) - return f, Py(r).to_numpy() -end - -mse(x::Matrix{T}, y::Matrix{T}) where {T<:Number} = _mse(x, y) - """ studentst(x, y) @@ -37,18 +30,11 @@ and its derivative w.r.t x `(k + 1) * (x - y) / (k + (x - y)^2)` """ -function _studentst(x::AT, y::AbstractArray{T}; k=T(2)) where {T<:Number, AT<:AbstractArray{T}} +function studentst(x::AbstractArray{T}, y::AbstractArray{T}; k=T(2)) where {T<:Number} k = convert(T, k) f = sum(_studentst_loss.(x, y, k)) r = (k + 1) .* (x - y) ./ (k .+ (x - y).^2) - return f::T, r::AT{T} -end - -function studentst(x::PyArray{T}, y::AbstractArray{T}; k=T(2)) where {T<:Number} - f, r = _studentst(x, y, k) - return f, Py(r).to_numpy() + return f, r end -studentst(x::Matrix{T}, y::Matrix{T}; k=T(2)) where {T<:Number} = _studentst(x, y, k) - _studentst_loss(x::T, y::T, k::T) where {T<:Number} = T(1/2) * (k + 1) * log(1 + (x-y)^2 / k) \ No newline at end of file diff --git a/src/TimeModeling/Modeling/misfit_fg.jl b/src/TimeModeling/Modeling/misfit_fg.jl index 3be8c0780..2eebfbc2c 100644 --- a/src/TimeModeling/Modeling/misfit_fg.jl +++ b/src/TimeModeling/Modeling/misfit_fg.jl @@ -65,9 +65,6 @@ function _multi_src_fg(model_full::AbstractModel, source::Dtypes, dObs::Dtypes, length(options.frequencies) == 0 ? freqs = nothing : freqs = options.frequencies - # Make dObserved/Y numpy to avoid indexing issues - dIn = Py(dObserved).to_numpy() - @juditime "Python call to J_adjoint" begin argout = wrapcall_data(ac.J_adjoint, modelPy, src_coords, qIn, rec_coords, dObserved, t_sub=options.subsampling_factor, checkpointing=options.optimal_checkpointing, diff --git a/src/TimeModeling/Modeling/python_interface.jl b/src/TimeModeling/Modeling/python_interface.jl index 1d5f69243..b518692e9 100644 --- a/src/TimeModeling/Modeling/python_interface.jl +++ b/src/TimeModeling/Modeling/python_interface.jl @@ -66,7 +66,6 @@ function devito_interface(modelPy::Py, srcGeometry::Nothing, srcData::Array, rec rec_coords = setup_grid(recGeometry, modelPy.shape) # Devito call - srcData = Py(srcData).to_numpy() return wrapcall_data(ac.forward_wf_src, modelPy, srcData, rec_coords, fw=fw, f0=options.f0, illum=illum) end @@ -77,7 +76,6 @@ function devito_interface(modelPy::Py, srcGeometry::Nothing, srcData::Array, rec dtComp = pyconvert(Float32, modelPy.critical_dt) # Devito call - srcData = Py(srcData).to_numpy() return wrapcall_data(ac.forward_wf_src_norec, modelPy, srcData, fw=fw, illum=illum) end @@ -115,9 +113,6 @@ function devito_interface(modelPy::Py, srcGeometry::Geometry, srcData::Array, re rec_coords = setup_grid(recGeometry, modelPy.shape) length(options.frequencies) == 0 ? freqs = nothing : freqs = options.frequencies - # Make dIn numpy to avoid indexing issues - dIn = Py(dIn).to_numpy() - return wrapcall_data(ac.J_adjoint, modelPy, src_coords, qIn, rec_coords, dIn, fw=fw, t_sub=options.subsampling_factor, checkpointing=options.optimal_checkpointing, @@ -131,16 +126,16 @@ end function devito_interface(modelPy::Py, weights::Array, srcData::Array, recGeometry::Geometry, recData::Nothing, dm::Nothing, options::JUDIOptions, illum::Bool, fw::Bool) judilog("Pr*$(_op_str(fw))*Pw'*w") - weights = pad_array(reshape(weights, modelPy.shape), modelPy.padsizes; mode=:zeros) + shape = pyconvert(Tuple, modelPy.shape) + weights = pad_array(reshape(weights, shape), modelPy.padsizes; mode=:zeros) # Interpolate input data to computational grid dtComp = pyconvert(Float32, modelPy.critical_dt) qIn = time_resample(srcData, recGeometry, dtComp) # Set up coordinates with devito dimensions - rec_coords = setup_grid(recGeometry, modelPy.shape) + rec_coords = setup_grid(recGeometry, shape) # Devito call - weights = Py(weights).to_numpy() return wrapcall_data(ac.forward_rec_w, modelPy, weights, qIn, rec_coords, fw=fw, f0=options.f0, illum=illum) end @@ -166,13 +161,14 @@ end function devito_interface(modelPy::Py, weights::Array, srcData::Array, recGeometry::Geometry, recData::Nothing, dm::Union{PhysicalParameter, Array}, options::JUDIOptions, illum::Bool, fw::Bool) judilog("Jw($(_op_str(fw)), q)*dm") - weights = pad_array(reshape(weights, modelPy.shape), modelPy.padsizes; mode=:zeros) + shape = pyconvert(Tuple, modelPy.shape) + weights = pad_array(reshape(weights, shape), modelPy.padsizes; mode=:zeros) # Interpolate input data to computational grid dtComp = pyconvert(Float32, modelPy.critical_dt) qIn = time_resample(srcData, recGeometry, dtComp) # Set up coordinates with devito dimensions - rec_coords = setup_grid(recGeometry, modelPy.shape) + rec_coords = setup_grid(recGeometry, shape) # Devito call return wrapcall_data(ac.born_rec_w, modelPy, weights, qIn, rec_coords, @@ -182,7 +178,8 @@ end # Adjoint Jacobian of extended source modeling: dm = J'*d_lin function devito_interface(modelPy::Py, weights::Array, srcData::Array, recGeometry::Geometry, recData::Array, dm::Nothing, options::JUDIOptions, illum::Bool, fw::Bool) judilog("Jw($(_op_str(fw)), q)'*d_lin") - weights = pad_array(reshape(weights, modelPy.shape), modelPy.padsizes; mode=:zeros) + shape = pyconvert(Tuple, modelPy.shape) + weights = pad_array(reshape(weights, shape), modelPy.padsizes; mode=:zeros) # Interpolate input data to computational grid dtComp = pyconvert(Float32, modelPy.critical_dt) qIn = time_resample(srcData, recGeometry, dtComp) @@ -190,12 +187,9 @@ function devito_interface(modelPy::Py, weights::Array, srcData::Array, recGeomet qIn, dIn = _maybe_pad_t0(qIn, recGeometry, dIn, recGeometry, dtComp) # Set up coordinates with devito dimensions - rec_coords = setup_grid(recGeometry, modelPy.shape) + rec_coords = setup_grid(recGeometry, shape) length(options.frequencies) == 0 ? freqs = nothing : freqs = options.frequencies - # Make dIn numpy to avoid indexing issues - dIn = Py(dIn).to_numpy() - return wrapcall_data(ac.J_adjoint, modelPy, nothing, qIn, rec_coords, dIn, fw=fw, t_sub=options.subsampling_factor, checkpointing=options.optimal_checkpointing, diff --git a/src/TimeModeling/Modeling/twri_objective.jl b/src/TimeModeling/Modeling/twri_objective.jl index 0f2ad059b..f64f2dd82 100644 --- a/src/TimeModeling/Modeling/twri_objective.jl +++ b/src/TimeModeling/Modeling/twri_objective.jl @@ -94,13 +94,6 @@ function _twri_objective(model_full::AbstractModel, source::judiVector, dObs::ju ~isempty(options.frequencies) ? freqs = options.frequencies : freqs = nothing ~isempty(options.frequencies) ? (wfilt, freqs) = filter_w(qIn, dtComp, freqs) : wfilt = nothing - # Make dObserved/Y numpy to avoid indexing issues - dObserved = Py(dObserved).to_numpy() - qIn = Py(qIn).to_numpy() - if !isnothing(Y) - Y = Py(Y).to_numpy() - end - argout = wrapcall_data(ac.wri_func, modelPy, src_coords, qIn, rec_coords, dObserved, Y, t_sub=options.subsampling_factor, grad=optionswri.params, grad_corr=optionswri.grad_corr, eps=optionswri.eps, diff --git a/src/TimeModeling/Preconditioners/DataPreconditioners.jl b/src/TimeModeling/Preconditioners/DataPreconditioners.jl index 59346a280..3b8a14ce6 100644 --- a/src/TimeModeling/Preconditioners/DataPreconditioners.jl +++ b/src/TimeModeling/Preconditioners/DataPreconditioners.jl @@ -118,7 +118,7 @@ function muteshot!(shot::AbstractMatrix{T}, rGeom::Geometry, srcGeom::Geometry; end end -function muteshot(shot::VecOrMat{T}, srcGeom::Geometry, recGeom::Geometry; +function muteshot(shot::AbstractVecOrMat{T}, srcGeom::Geometry, recGeom::Geometry; vp=1500, t0=.1, mode=:reflection, taperwidth=floor(Int, 2/t0)) where {T<:Number} sr = reshape(shot, recGeom) for s=1:get_nsrc(recGeom) @@ -158,7 +158,7 @@ judiFilter(geometry::Geometry, fmin::T, fmax::T) where T = judiFilter(geometry, judiFilter(geometry::Geometry, fmin::Float32, fmax::Float32) = FrequencyFilter{Float32, fmin, fmax}(n_samples(geometry), geometry) judiFilter(v::judiVector, fmin, fmax) = judiFilter(v.geometry, fmin, fmax) -function matvec(D::FrequencyFilter{T, fm, FM}, x::VecOrMat{T}) where {T, fm, FM} +function matvec(D::FrequencyFilter{T, fm, FM}, x::AbstractVecOrMat{T}) where {T, fm, FM} dr = reshape(x, D.recGeom) for j=1:get_nsrc(D.recGeom) dri = view(dr, :, :, j) @@ -227,7 +227,7 @@ filter_data(Din::judiVector, ::Any; fmin=0, fmax=Inf) = filter_data(Din; fmin=fm Performs a causal filtering [fmin, fmax] on the input data bases on its sampling rate `dt`. Automatically perfroms a lowpass if fmin=0 (default) """ -function filter_data(Din::Matrix{T}, dt_in; fmin=0, fmax=Inf) where T +function filter_data(Din::AbstractMatrix{T}, dt_in; fmin=0, fmax=Inf) where T out = similar(Din) filter!(out, Din, dt_in; fmin=T(fmin), fmax=T(fmax)) return out @@ -304,7 +304,7 @@ function matvec(D::TimeDifferential{T, K}, x::judiVector{T, AT}) where {T, AT, K return out end -function matvec(D::TimeDifferential{T, K}, x::Array{T}) where {T, K} +function matvec(D::TimeDifferential{T, K}, x::AbstractArray{T}) where {T, K} xr = reshape(x, D.recGeom) out = similar(xr) # make omega^K @@ -371,7 +371,7 @@ function matvec(D::TimeGain{T, K}, x::judiVector{T, AT}) where {T, AT, K} return out end -function matvec(D::TimeGain{T, K}, x::Array{T}) where {T, K} +function matvec(D::TimeGain{T, K}, x::AbstractArray{T}) where {T, K} xr = reshape(x, D.recGeom) # make time^K timek = (D.recGeom.taxis[1]).^K diff --git a/src/TimeModeling/Preconditioners/base.jl b/src/TimeModeling/Preconditioners/base.jl index 8662f6f22..94ebfb96d 100644 --- a/src/TimeModeling/Preconditioners/base.jl +++ b/src/TimeModeling/Preconditioners/base.jl @@ -28,6 +28,7 @@ getproperty(J::Preconditioner, s::Symbol) = _get_property(J, Val{s}()) *(J::Preconditioner, ms::judiMultiSourceVector) = matvec(J, ms) *(J::Preconditioner, ms::PhysicalParameter) = matvec(J, ms) *(J::Preconditioner, v::VecOrMat{T}) where T = matvec(J, v) +*(J::Preconditioner, v::PyArray{T}) where T = matvec(J, v) mul!(out::judiMultiSourceVector, J::Preconditioner, ms::judiMultiSourceVector) = copyto!(out, matvec(J, ms)) mul!(out::PhysicalParameter, J::Preconditioner, ms::PhysicalParameter) = copyto!(out, matvec(J, ms)) diff --git a/src/TimeModeling/Utils/auxiliaryFunctions.jl b/src/TimeModeling/Utils/auxiliaryFunctions.jl index e40b6b41e..ab3910c37 100644 --- a/src/TimeModeling/Utils/auxiliaryFunctions.jl +++ b/src/TimeModeling/Utils/auxiliaryFunctions.jl @@ -25,7 +25,8 @@ Parameters """ function devito_model(model::MT, options::JUDIOptions, dm) where {MT<:AbstractModel} # Set up Python model structure - physpar = Dict((n, isa(v, PhysicalParameter) ? Py(v.data).to_numpy() : v) for (n, v) in _params(model)) + physpar = Dict((n, isa(v, PhysicalParameter) ? v.data : v) for (n, v) in _params(model)) + dm = isnothing(dm) ? nothing : (isa(dm, PhysicalParameter) ? dm.data : dm) modelPy = pm.Model(origin(model), spacing(model), size(model), fs=options.free_surface, nbl=nbl(model), space_order=options.space_order, dt=options.dt_comp, dm=dm; @@ -52,16 +53,17 @@ function pad_array(m::Array{DT}, nb::NTuple{N, Tuple{Int64, Int64}}; mode::Symbo n = size(m) new_size = Tuple([n[i] + sum(nb[i]) for i=1:length(nb)]) Ei = [] - for i=1:length(nb) + for i=length(nb):-1:1 left, right = nb[i] push!(Ei, joExtend(n[i], mode;pad_upper=right, pad_lower=left, RDT=DT, DDT=DT)) end - padded = joKron(Ei...) * PermutedDimsArray(m, length(n):-1:1)[:] - return PyReverseDims(reshape(padded, reverse(new_size))) + padded = joKron(Ei...) * m[:] + return reshape(padded, new_size) end pad_array(::Nothing, ::NTuple{N, Tuple{Int64, Int64}}; s::Symbol=:border) where N = nothing pad_array(m::Number, ::NTuple{N, Tuple{Int64, Int64}}; s::Symbol=:border) where N = m +pad_array(m::Array{DT}, nb::Py; mode::Symbol=:border) where DT = pad_array(m, pyconvert(Tuple, nb), mode=mode) """ remove_padding(m, nb; true_adjoint=False) @@ -86,7 +88,7 @@ function remove_padding(gradient::AbstractArray{DT}, nb::NTuple{ND, Tuple{Int64, return out end -remove_padding(g::PyArray, nb::Py; true_adjoint::Bool=false) = remove_padding(g, pyconvert(Tuple, nb), true_adjoint=true) +remove_padding(g::PyArray, nb::Py; true_adjoint::Bool=false) = remove_padding(g, pyconvert(Tuple, nb), true_adjoint=true_adjoint) """ limit_model_to_receiver_area(srcGeometry, recGeometry, model, buffer; pert=nothing) @@ -390,12 +392,6 @@ setup_grid(geometry::GeometryIC{T}, ::NTuple{2, <:Integer}) where {T<:Real} = hc setup_grid(geometry::GeometryIC{T}, ::NTuple{1, <:Integer}) where {T<:Real} = geometry.xloc[1] setup_grid(geometry::GeometryIC{T}, t::Py) where T<:Real = setup_grid(geometry, pyconvert(Tuple, t)) -""" - setup_3D_grid(x, y, z) - -""" - setup_3D_grid(x, y, z) - """ setup_3D_grid(x, y, z) diff --git a/src/TimeModeling/Utils/time_utilities.jl b/src/TimeModeling/Utils/time_utilities.jl index 4deaa2812..8b008f59a 100644 --- a/src/TimeModeling/Utils/time_utilities.jl +++ b/src/TimeModeling/Utils/time_utilities.jl @@ -76,8 +76,10 @@ function time_resample(data::AbstractArray{T, N}, dt_in::Real, G_in::Geometry{T} return dout end +_time_resample(data::PyArray{T, 2}, rate::Integer) where T = data[1:rate:end, :] _time_resample(data::Matrix{T}, rate::Integer) where T = data[1:rate:end, :] _time_resample(data::PermutedDimsArray{T, 2, (2, 1), (2, 1), Matrix{T}}, rate::Integer) where {T<:Real} = data.parent[:, 1:rate:end]' +_time_resample(data::PermutedDimsArray{T, 2, (2, 1), (2, 1), PyArray{T, 2}}, rate::Integer) where {T<:Real} = data.parent[:, 1:rate:end]' SincInterpolation(Y::AbstractMatrix{T}, S::StepRangeLen{T}, Up::StepRangeLen{T}) where T<:Real = sinc.( (Up .- S') ./ (S[2] - S[1]) ) * Y diff --git a/src/pysource/fields.py b/src/pysource/fields.py index 98f5c8870..66ef30f4e 100644 --- a/src/pysource/fields.py +++ b/src/pysource/fields.py @@ -100,7 +100,7 @@ def src_wavefield(model, u, fw=True): Forward or backward (for naming) """ name = "uqwf" if fw else "vqwf" - init = u.data if isinstance(u, TimeFunction) else u + init = u.data if isinstance(u, TimeFunction) else u.to_numpy(copy=False) wf_src = TimeFunction(name=name, grid=model.grid, time_order=2, space_order=0, save=u.shape[0], initializer=init) return wf_src @@ -179,11 +179,12 @@ def lr_src_fields(model, weight, wavelet, empty_w=False, rec=False): wavelett = TimeFunction(name='wf_%s' % wn, dimensions=(time,), time_dim=time, shape=(nt,), save=nt, grid=model.grid) wavelett.data[:] = np.array(wavelet)[:, 0] - if empty_w: - source_weight = Function(name='%s_weight' % wn, grid=model.grid, space_order=0) + if isinstance(weight, Function): + source_weight = weight else: source_weight = Function(name='%s_weight' % wn, grid=model.grid, space_order=0) - initialize_function(source_weight, weight, 0) + if not empty_w: + initialize_function(source_weight, weight.to_numpy(copy=False), 0) return source_weight, wavelett diff --git a/src/pysource/geom_utils.py b/src/pysource/geom_utils.py index 522da3b5e..6ba0ee740 100644 --- a/src/pysource/geom_utils.py +++ b/src/pysource/geom_utils.py @@ -27,7 +27,8 @@ def src_rec(model, u, src_coords=None, rec_coords=None, wavelet=None, nt=None): return src, rcv -def geom_expr(model, u, src_coords=None, rec_coords=None, wavelet=None, fw=True, nt=None): +def geom_expr(model, u, src_coords=None, rec_coords=None, + wavelet=None, fw=True, nt=None): """ Generates the source injection and receiver interpolation. This function is fully abstracted and does not care whether this is a diff --git a/src/pysource/interface.py b/src/pysource/interface.py index 15513bc9f..4ed2c69c1 100644 --- a/src/pysource/interface.py +++ b/src/pysource/interface.py @@ -8,7 +8,7 @@ from propagators import forward, born, gradient, forward_grad from sensitivity import Loss from sources import Receiver -from utils import weight_fun, compute_optalpha +from utils import weight_fun, compute_optalpha, npdot from fields import memory_field, src_wavefield from fields_exprs import wf_as_src @@ -584,7 +584,7 @@ def wri_func(model, src_coords, wavelet, rec_coords, recin, yin, if yin is None or grad_corr: y, u0, _, _ = forward(model, src_coords, rec_coords, wavelet, save=grad_corr, ws=ws, f0=f0) - ydat = recin[:] - y.data[:] + ydat = recin - y.data[:] else: ydat = yin @@ -596,8 +596,7 @@ def wri_func(model, src_coords, wavelet, rec_coords, recin, yin, c2 = np.log(np.prod(model.shape)) # = - ndt = np.sqrt(model.critical_dt) - PTy_dot_r = ndt**2 * (np.dot(ydat.reshape(-1), recin.reshape(-1)) - - np.dot(srca.data.reshape(-1), wavelet.reshape(-1))) + PTy_dot_r = ndt**2 * (npdot(ydat, recin) - npdot(srca.data, wavelet)) norm_y = ndt * np.linalg.norm(ydat) # alpha diff --git a/src/pysource/models.py b/src/pysource/models.py index b811d5927..da2965909 100644 --- a/src/pysource/models.py +++ b/src/pysource/models.py @@ -41,6 +41,13 @@ def getmax(f): return np.max(f) +def to_numpy(f): + try: + return f.to_numpy(copy=False) + except AttributeError: + return f + + _thomsen = [('epsilon', 1), ('delta', 1), ('theta', 0), ('phi', 0)] @@ -208,8 +215,8 @@ def __init__(self, origin, spacing, shape, space_order=8, nbl=40, dtype=np.float # Additional parameter fields for TTI operators if self._is_tti: - epsilon = 1 if epsilon is None else 1 + 2 * epsilon - delta = 1 if delta is None else 1 + 2 * delta + epsilon = 1 if epsilon is None else 1 + 2 * to_numpy(epsilon) + delta = 1 if delta is None else 1 + 2 * to_numpy(delta) self.epsilon = self._gen_phys_param(epsilon, 'epsilon', space_order) self.scale = np.sqrt(np.max(epsilon)) self.delta = self._gen_phys_param(delta, 'delta', space_order) @@ -221,8 +228,9 @@ def __init__(self, origin, spacing, shape, space_order=8, nbl=40, dtype=np.float if self._is_elastic: self.lam = self._gen_phys_param(lam, 'lam', space_order, is_param=True) try: + mu = to_numpy(mu) mu[np.where(mu == 0)] = 1e-12 - except TypeError: + except (ValueError, TypeError): mu = 1e-12 if mu == 0 else mu self.mu = self._gen_phys_param(mu, 'mu', space_order, is_param=True, avg_mode='harmonic') @@ -230,7 +238,7 @@ def __init__(self, origin, spacing, shape, space_order=8, nbl=40, dtype=np.float b = b if b is not None else 1 / rho except TypeError: b = 1 - vp_vals = ((lam + 2 * mu) * b)**(.5) + vp_vals = ((to_numpy(lam) + 2 * mu) * b)**(.5) # User provided dt self._dt = kwargs.get('dt') @@ -297,8 +305,9 @@ def _gen_phys_param(self, field, name, space_order, is_param=False, """ if field is None: return default_value - if isinstance(field, np.ndarray): - if _dtypes['params'] == 'f16' or field.dtype == np.float16: + if hasattr(field, '__array_interface__') and field.shape: + dtype = np.dtype(field.__array_interface__['typestr']).type + if _dtypes['params'] == 'f16' or dtype == np.float16: _min = np.amin(field) _max = np.amax(field) if _max == _min: @@ -310,7 +319,7 @@ def _gen_phys_param(self, field, name, space_order, is_param=False, function = Function(name=name, grid=self.grid, space_order=space_order, parameter=is_param, avg_mode=avg_mode, dtype=dtype) pad = 0 if field.shape == function.shape else self.padsizes - initialize_function(function, field, pad) + initialize_function(function, to_numpy(field), pad) else: function = Constant(name=name, value=np.amin(field)) self._physical_parameters.append(name) @@ -497,7 +506,7 @@ def abox(self, src, rec, fw=True): abox = ABox(src, rec, self._vp, self.space_order, eps=eps)._abox_func abox.data[:] = abox._compute(src=src, rcv=rec, vp=self._vp, eps=eps, dt=self.critical_dt) - return {'abox': abox} + return {'srcbox': abox} def __init_abox__(self, src, rec, fw=True): return @@ -586,4 +595,4 @@ def __init_abox__(self, src, rec, fw=True): if not fw: src, rec = rec, src eps = getattr(self, 'epsilon', None) - self._abox = ABox(src, rec, self._vp, self.space_order, eps=eps) + self._abox = ABox(src, rec, self._vp, self.space_order, eps=eps, name='srcbox') diff --git a/src/pysource/propagators.py b/src/pysource/propagators.py index c6ebc1978..97a300872 100644 --- a/src/pysource/propagators.py +++ b/src/pysource/propagators.py @@ -214,7 +214,7 @@ def born(model, src_coords, rcv_coords, wavelet, save=False, # Update kwargs kw.update(fields_kwargs(u, ul, snl, rnl, rcvl, u_save, dft_m, fr, ws, wt, f0q, I)) kw.update(model.physical_params(born=True)) - kw.update(model.abox(snl, rnl, fw=True)) + kw.update(model.abox(snl, rcvl, fw=True)) # SLS field if model.is_viscoacoustic: diff --git a/src/pysource/sensitivity.py b/src/pysource/sensitivity.py index 1558fe16d..9f5ca0326 100644 --- a/src/pysource/sensitivity.py +++ b/src/pysource/sensitivity.py @@ -19,7 +19,7 @@ def func_name(freq=None, ic="as"): Get key for imaging condition/linearized source function """ if freq is None: - return ic + return str(ic) else: return "%s_%s" % (ic, "freq") @@ -170,7 +170,7 @@ def lin_src(model, u, ic="as"): ic: String Imaging condition of which we compute the linearized source """ - ls_func = ls_dict[func_name(ic=ic)] + ls_func = ls_dict[func_name(ic=str(ic))] return ls_func(model, as_tuple(u)) @@ -258,11 +258,11 @@ def Loss(dsyn, dobs, dt, is_residual=False, misfit=None): if misfit is not None: if isinstance(dsyn, tuple): f, r = misfit(dsyn[0].data._local, dobs - dsyn[1].data._local[:]) - dsyn[0].data._local[:] = r[:] + dsyn[0].data._local[:] = r return dt * f, dsyn[0].data._local else: f, r = misfit(dsyn.data._local, dobs) - dsyn.data._local[:] = r[:] + dsyn.data._local[:] = r return dt * f, dsyn.data._local if not is_residual: diff --git a/src/pysource/utils.py b/src/pysource/utils.py index 01ddf39be..4809879e7 100644 --- a/src/pysource/utils.py +++ b/src/pysource/utils.py @@ -160,3 +160,11 @@ def compression_mode(): return 'bitcomp' else: return 'noop' + + +def npdot(a, b): + """ + Inner product of n-dimensional ndarrays through numpy einsum + """ + inds = 'ijklm'[:len(a.shape)] + return np.einsum('%s,%s->' % (inds, inds), a, b) diff --git a/src/utilities.jl b/src/utilities.jl index 0b3b5be7b..55231a6dd 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -59,8 +59,8 @@ end # Devito configuration -set_devito_config(key::String, val::String) = devito[].configuration[key] = val -set_devito_config(key::String, val::Bool) = devito[].configuration[key] = val +set_devito_config(key::String, val::String) = begin devito.configuration[key] = val end +set_devito_config(key::String, val::Bool) = begin devito.configuration[key] = val end set_devito_config(kw...) = begin for (k, v) in kw diff --git a/test/test_basics.jl b/test/test_basics.jl index 9195b4483..7d07a9e95 100644 --- a/test/test_basics.jl +++ b/test/test_basics.jl @@ -51,13 +51,13 @@ function test_density(ndim) model = Model(n, d, o, m, rho; nb=0) @test :rho in JUDI._mparams(model) modelpy = devito_model(model, Options()) - @test isapprox(modelpy.irho.data, 1 ./ model.rho) + @test isapprox(JUDI.PythonCall.PyArray(modelpy.irho.data), 1 ./ model.rho) rho[61] = 1000 model = Model(n, d, o, m, rho; nb=0) @test :rho in JUDI._mparams(model) modelpy = devito_model(model, Options()) - @test isapprox(modelpy.rho.data, model.rho) + @test isapprox(JUDI.PythonCall.PyArray(modelpy.rho.data), model.rho) end @@ -68,11 +68,13 @@ function test_padding(ndim) modelPy = devito_model(model, Options()) m0 = model.m - mcut = remove_padding(deepcopy(modelPy.m.data), modelPy.padsizes; true_adjoint=true) - mdata = deepcopy(modelPy.m.data) + mdata = Array(JUDI.PythonCall.PyArray(modelPy.m.data)) + mcut = remove_padding(copy(mdata), JUDI.PythonCall.pyconvert(Tuple, modelPy.padsizes); true_adjoint=true) - @show dot(m0, mcut)/dot(mdata, mdata) - @test isapprox(dot(m0, mcut), dot(mdata, mdata)) + a = dot(m0, mcut) + b = dot(mdata, mdata) + @show a/b, (a - b) / (a + b) + @test (a - b) / (a + b) ≈ 0 atol=5e-6 rtol=0 end function test_limit_m(ndim, tti) @@ -88,8 +90,8 @@ function test_limit_m(ndim, tti) srcGeometry = example_src_geometry() recGeometry = example_rec_geometry(cut=true) buffer = 100f0 - dm = rand(Float32, model.n...) - dmb = PhysicalParameter(0 .* dm, model.n, model.d, model.o) + dm = rand(Float32, size(model)...) + dmb = PhysicalParameter(0 .* dm, size(model), spacing(model), origin(model)) new_mod, dm_n = limit_model_to_receiver_area(srcGeometry, recGeometry, deepcopy(model), buffer; pert=dm) # check new_mod coordinates @@ -120,9 +122,9 @@ function test_limit_m(ndim, tti) dmt = PhysicalParameter(dm_n, new_mod.n, model.d, new_mod.o) ex_dm = dmb + dmt - @test ex_dm.o == model.o - @test ex_dm.n == model.n - @test ex_dm.d == model.d + @test ex_dm.o == origin(model) + @test ex_dm.n == size(model) + @test ex_dm.d == spacing(model) @test norm(ex_dm) == norm(dm_n) end @@ -177,14 +179,24 @@ function test_serial() end @testset "Test basic utilities" begin - @timeit TIMEROUTPUT "Basic setup utilities" begin - test_ftp() + # @timeit TIMEROUTPUT "FTP data" begin + # test_ftp() + # end + @timeit TIMEROUTPUT "Serial/parallel" begin test_serial() + end + @timeit TIMEROUTPUT "3D setup" begin setup_3d() + end + + @timeit TIMEROUTPUT "Padding and density" begin for ndim=[2, 3] test_padding(ndim) test_density(ndim) end + end + + @timeit TIMEROUTPUT "Limit m" begin opt = Options(frequencies=[[2.5, 4.5], [3.5, 5.5]]) @test subsample(opt, 1).frequencies == [2.5, 4.5] @test subsample(opt, 2).frequencies == [3.5, 5.5] @@ -194,7 +206,9 @@ end test_limit_m(ndim, tti) end end + end + @timeit TIMEROUTPUT "Model config" begin # Test model for ndim=[2, 3] for (tti, elas, visco) in [(false, false, false), (true, false, false), (false, true, false), (false, false, true)] @@ -202,13 +216,13 @@ end # Default dt modelPy = devito_model(model, Options()) - @test isapprox(modelPy.critical_dt, calculate_dt(model)) + @test isapprox(JUDI.PythonCall.pyconvert(Float32, modelPy.critical_dt), calculate_dt(model)) @test get_dt(model) == calculate_dt(model) @test isapprox(calculate_dt(model; dt=.5f0), .5f0) # Input dt modelPy = devito_model(model, Options(dt_comp=.5f0)) - @test modelPy.critical_dt == .5f0 + @test Bool(modelPy.critical_dt == .5f0) # Verify nt srcGeometry = example_src_geometry() @@ -225,4 +239,5 @@ end end end end + end \ No newline at end of file diff --git a/test/test_preconditioners.jl b/test/test_preconditioners.jl index bbf0bc9bb..21aac05c4 100644 --- a/test/test_preconditioners.jl +++ b/test/test_preconditioners.jl @@ -3,7 +3,7 @@ model, model0, dm = setup_model(tti, viscoacoustic, nlayer) wb = find_water_bottom(model.m .- maximum(model.m)) -q, srcGeometry, recGeometry = setup_geom(model; nsrc=2) +q, srcGeometry, recGeometry = setup_geom(model; nsrc=2, tn=1750f0) ftol = 1f-5 @@ -203,7 +203,7 @@ dm = model0.m - model.m # Illumination I = judiIllumination(FM; mode="v") - dloc = FM*q + dloc = FM[1]*q[1] # forward only, nothing done @test I.illums["v"].data == ones(Float32, model.n) @@ -219,6 +219,8 @@ dm = model0.m - model.m @test "v" ∈ keys(Iv.illums) @test "u" ∉ keys(Iv.illums) @test norm(Iv.illums["v"]) == norm(ones(Float32, model.n)) + @test I.illums["u"].data == bck + # Test Product @test inv(I)*I*model0.m ≈ model0.m rtol=ftol atol=0