Merge pull request #157 from slimgroup/issue-156

Fixes padding for limit_m

.github/workflows/ci-judi.yml

@@ -50,6 +50,11 @@ jobs:
         if: runner.os == 'macOS'
         run: brew install gcc@9
 
+      - name: Set up Python 3.9
+        uses: actions/setup-python@v4
+        with:
+          python-version: 3.9
+
       - name: Set julia python
         run: |
           PYTHON=$(which python3) julia -e 'using Pkg;Pkg.add("PyCall");Pkg.build("PyCall")'

.github/workflows/ci-op.yml

@@ -79,6 +79,11 @@ jobs:
           version: ${{ matrix.version }}
           arch: x64
 
+      - name: Set up Python 3.9
+        uses: actions/setup-python@v4
+        with:
+          python-version: 3.9
+
       - name: Set julia python
         run: |
           PYTHON=$(which python3) julia -e 'using Pkg;Pkg.add("PyCall");Pkg.build("PyCall")'

.github/workflows/deploy_doc.yaml

@@ -16,7 +16,12 @@ jobs:
     steps:
       - name: Checkout master
         uses: actions/checkout@v2
-
+
+      - name: Set up Python 3.9
+        uses: actions/setup-python@v4
+        with:
+          python-version: 3.9
+
       - uses: julia-actions/setup-julia@latest
 
       - name: Install dependencies

Project.toml

@@ -1,7 +1,7 @@
 name = "JUDI"
 uuid = "f3b833dc-6b2e-5b9c-b940-873ed6319979"
 authors = ["Philipp Witte, Mathias Louboutin"]
-version = "3.1.12"
+version = "3.1.13"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

src/TimeModeling/Modeling/misfit_fg.jl

@@ -108,7 +108,7 @@ Example
 """
 function fwi_objective(model::MTypes, q::Dtypes, dobs::Dtypes; options=Options(), kw...)
     n_exp = check_args(model, q, dobs)
-    G = n_exp == 1 ? similar(model.m) : [similar(get_exp(model, i).m) for i=1:n_exp]
+    G = n_exp == 1 ? similar(model.m, model) : [similar(get_exp(model, i).m, get_exp(model, i)) for i=1:n_exp]
     f = fwi_objective!(G, model, q, dobs; options=options, kw...)
     f, G
 end
@@ -126,7 +126,7 @@ Example
 """
 function lsrtm_objective(model::MTypes, q::Dtypes, dobs::Dtypes, dm::dmTypes; options=Options(), nlind=false, kw...)
     n_exp = check_args(model, q, dobs, dm)
-    G = n_exp == 1 ? similar(model.m) : [similar(get_exp(model, i).m) for i=1:n_exp]
+    G = n_exp == 1 ? similar(model.m, model) : [similar(get_exp(model, i).m, get_exp(model, i)) for i=1:n_exp]
     f = lsrtm_objective!(G, model, q, dobs, dm; options=options, nlind=nlind, kw...)
     f, G
 end

src/TimeModeling/Modeling/propagation.jl

@@ -66,6 +66,7 @@ function multi_src_propagate(F::judiPropagator{T, O}, q::AbstractVector{T}) wher
     arg_func = i -> (F[i], src_i(F, q, i))
     # Distribute source
     res = run_and_reduce(propagate, pool, nsrc, arg_func)
+    res = _project_to_physical_domain(res, F.model)
     res = as_vec(res, Val(F.options.return_array))
     return res
 end
@@ -87,7 +88,7 @@ function multi_src_fg!(G, model, q, dobs, dm; options=Options(), kw...)
     # Distribute source
     res = run_and_reduce(multi_src_fg, pool, nsrc, arg_func)
     f, g = as_vec(res, Val(options.return_array))
-    G .= g
+    G .+= g
     return f
 end
 

src/TimeModeling/Modeling/time_modeling_serial.jl

@@ -13,12 +13,8 @@ function time_modeling(model_full::Model, srcGeometry::GeomOrNot, srcData::Array
     srcGeometry = Geometry(srcGeometry)
 
     # Reutrn directly for J*0
-    if op==:born
-        if norm(dm) == 0 && options.return_array == false
-            return judiVector(recGeometry, zeros(Float32, recGeometry.nt[1], length(recGeometry.xloc[1])))
-        elseif norm(dm) == 0 && options.return_array == true
-            return vec(zeros(Float32, recGeometry.nt[1], length(recGeometry.xloc[1])))
-        end
+    if (op==:born && norm(dm) == 0)
+        return judiVector(recGeometry, zeros(Float32, recGeometry.nt[1], length(recGeometry.xloc[1])))
     end
 
     # limit model to area with sources/receivers

src/TimeModeling/Types/ModelStructure.jl

@@ -340,4 +340,14 @@ ndims(m::Model) = ndims(m.m.data)
 
 display(m::Model) = println("Model (n=$(m.n), d=$(m.d), o=$(m.o)) with parameters $(keys(m.params))")
 show(io::IO, m::Model) = print(io, "Model (n=$(m.n), d=$(m.d), o=$(m.o)) with parameters $(keys(m.params))")
-show(io::IO, ::MIME{Symbol("text/plain")}, m::Model) = print(io, "Model (n=$(m.n), d=$(m.d), o=$(m.o)) with parameters $(keys(m.params))")
+show(io::IO, ::MIME{Symbol("text/plain")}, m::Model) = print(io, "Model (n=$(m.n), d=$(m.d), o=$(m.o)) with parameters $(keys(m.params))")
+
+# Pad gradient if aperture doesn't match full domain
+_project_to_physical_domain(p, ::Any) = p
+
+function _project_to_physical_domain(p::PhysicalParameter, model::Model)
+    p.n == model.n && (return p)
+    pp = similar(p, model)
+    pp .+= p
+    return pp
+end

test/test_issues.jl

@@ -70,4 +70,51 @@ end
         dobs = F*q
         @test ~isnan(norm(dobs))
     end
+end
+
+@testset "Tests limit_m issue 156" begin
+    @timeit TIMEROUTPUT "Issue 156" begin
+        model, model0, dm = setup_model(tti, viscoacoustic, nlayer)
+        q, srcGeometry, recGeometry, f0 = setup_geom(model; nsrc=1)
+        # Restrict rec to middle of the model
+        recGeometry.xloc[1] .= range(11*model.d[1], stop=(model.n[1] - 12)*model.d[1],
+                                     length=length(recGeometry.xloc[1]))
+        # Data
+        F = judiModeling(model, srcGeometry, recGeometry)
+        dobs = F*q
+        # Run gradient and check output size
+        opt = Options(limit_m=true, buffer_size=0f0)
+        F0 = judiModeling(model0, srcGeometry, recGeometry; options=opt)
+
+        # fwi wrapper
+        g_ap = JUDI.multi_src_fg(model0, q , dobs, nothing, opt, false, false, mse)[2]
+        @test g_ap.n == (model.n .- (22, 0))
+        @test g_ap.o[1] == model.d[1]*11
+
+        g1 = fwi_objective(model0, q, dobs; options=opt)[2]
+        @test g1.n  == model.n
+        @test norm(g1.data[1:10, :]) == 0
+        @test norm(g1.data[end-10:end, :]) == 0
+
+        # lsrtm wrapper
+        g_ap = JUDI.multi_src_fg(model0, q , dobs, dm, opt, false, true, mse)[2]
+        @test g_ap.n == (model.n .- (22, 0))
+        @test g_ap.o[1] == model.d[1]*11
+
+        g2 = lsrtm_objective(model0, q, dobs, dm; options=opt)[2]
+        @test g2.n  == model.n
+        @test norm(g2.data[1:10, :]) == 0
+        @test norm(g2.data[end-10:end, :]) == 0
+
+        # Lin op
+        Jp = judiJacobian(F0, q)'
+        g_ap = JUDI.propagate(Jp, dobs)
+        @test g_ap.n == (model.n .- (22, 0))
+        @test g_ap.o[1] == model.d[1]*11
+
+        g3 = Jp * dobs
+        @test g3.n  == model.n
+        @test norm(g3.data[1:10, :]) == 0
+        @test norm(g3.data[end-10:end, :]) == 0
+    end
 end

test/test_jacobian.jl

@@ -8,31 +8,38 @@
 
 ### Model
 model, model0, dm = setup_model(tti, viscoacoustic, nlayer)
-q, srcGeometry, recGeometry, f0 = setup_geom(model)
+q, srcGeometry, recGeometry, f0 = setup_geom(model; nsrc=2)
 dt = srcGeometry.dt[1]
 
 m0 = model0.m
 ######################## WITH DENSITY ############################################
 
 @testset "Jacobian test with $(nlayer) layers and tti $(tti) and viscoacoustic $(viscoacoustic) freesurface $(fs)" begin
-    # Write shots as segy files to disk
-    opt = Options(sum_padding=true, dt_comp=dt, free_surface=fs, f0=f0)
+    @timeit TIMEROUTPUT "Jacobian generic tests" begin
+        # Write shots as segy files to disk
+        opt = Options(sum_padding=true, dt_comp=dt, free_surface=fs, f0=f0)
 
-    # Setup operators
-    F = judiModeling(model, srcGeometry, recGeometry; options=opt)
-    F0 = judiModeling(model0, srcGeometry, recGeometry; options=opt)
-    J = judiJacobian(F0, q)
-
-    # Linear modeling
-    dobs = F*q
-    dD = J*dm
+        # Setup operators
+        F = judiModeling(model, srcGeometry, recGeometry; options=opt)
+        F0 = judiModeling(model0, srcGeometry, recGeometry; options=opt)
+        J = judiJacobian(F0, q)
 
-    @test norm(J*(0f0.*dm)) == 0
-    @test isapprox(dD, J*vec(dm.data); rtol=1f-6)
+        # Linear modeling
+        dobs = F*q
+        dD = J*dm
+        dlin0 = J*(0f0.*dm)
+
+        @test norm(dlin0) == 0
+        @test isapprox(dD, J*vec(dm.data); rtol=1f-6)
 
-    # Gradient test
-    grad_test(x-> F(;m=x)*q, m0, dm, dD; data=true)
+        # Gradient test
+        grad_test(x-> F(;m=x)*q, m0, dm, dD; data=true)
 
+        # Check return_array returns correct size with zero dm and multiple shots
+        J.options.return_array = true
+        dlin0v = J*(0f0.*dm)
+        @test length(dlin0) == length(vec(dlin0))
+    end
 end
 
 ### Extended source
@@ -45,14 +52,14 @@ end
         Pr = judiProjection(recGeometry)
         F = judiModeling(model; options=opt)
         F0 = judiModeling(model0; options=opt)
-        Pw = judiLRWF(q.geometry.dt[1], q.data[1])
+        Pw = judiLRWF(q.geometry.dt, q.data)
 
         # Combined operators
         A = Pr*F*adjoint(Pw)
         A0 = Pr*F0*adjoint(Pw)
 
         # Extended source weights
-        w = judiWeights(randn(Float32, model0.n))
+        w = judiWeights(randn(Float32, model0.n); nsrc=2)
         J = judiJacobian(A0, w)
 
         # Nonlinear modeling