add combined J'J for muted data backpropagated input

Manifest.toml

@@ -18,17 +18,17 @@ git-tree-sha1 = "e928ca0a49f7b0564044b39108c70c160f03e05a"
 uuid = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 version = "1.1.2"
 
+[[ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+
 [[ArrayInterface]]
 deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
 git-tree-sha1 = "ce17bad65d0842b34a15fffc8879a9f68f08a67f"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 version = "3.1.6"
 
 [[Artifacts]]
-deps = ["Pkg"]
-git-tree-sha1 = "c30985d8821e0cd73870b17b0ed0ce6dc44cb744"
 uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
-version = "1.3.0"
 
 [[AxisAlgorithms]]
 deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
@@ -69,10 +69,10 @@ uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
 version = "0.4.1"
 
 [[CUDA]]
-deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CompilerSupportLibraries_jll", "DataStructures", "ExprTools", "GPUArrays", "GPUCompiler", "LLVM", "Libdl", "LinearAlgebra", "Logging", "MacroTools", "NNlib", "Pkg", "Printf", "Random", "Reexport", "Requires", "SparseArrays", "Statistics", "TimerOutputs"]
-git-tree-sha1 = "6ccc73b2d8b671f7a65c92b5f08f81422ebb7547"
+deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CompilerSupportLibraries_jll", "DataStructures", "ExprTools", "GPUArrays", "GPUCompiler", "LLVM", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "MacroTools", "Memoize", "NNlib", "Printf", "Random", "Reexport", "Requires", "SparseArrays", "Statistics", "TimerOutputs"]
+git-tree-sha1 = "870e029382294443a6578190e992bf4cbfd34e22"
 uuid = "052768ef-5323-5732-b1bb-66c8b64840ba"
-version = "2.4.1"
+version = "2.6.2"
 
 [[CatIndices]]
 deps = ["CustomUnitRanges", "OffsetArrays"]
@@ -82,9 +82,9 @@ version = "0.2.2"
 
 [[ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "1132f4e0893f286b1b88eb64ec5ddc8a5ae4b04e"
+git-tree-sha1 = "644c24cd6344348f1c645efab24b707088be526a"
 uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "0.9.32"
+version = "0.9.34"
 
 [[ColorTypes]]
 deps = ["FixedPointNumbers", "Random"]
@@ -117,10 +117,8 @@ uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
 version = "3.25.0"
 
 [[CompilerSupportLibraries_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "8e695f735fca77e9708e795eda62afdb869cbb70"
+deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
-version = "0.3.4+0"
 
 [[ComputationalResources]]
 git-tree-sha1 = "52cb3ec90e8a8bea0e62e275ba577ad0f74821f7"
@@ -209,6 +207,10 @@ git-tree-sha1 = "ab2f313a5ef4db9d127b8fae4d922e268db1552f"
 uuid = "aaf54ef3-cdf8-58ed-94cc-d582ad619b94"
 version = "0.6.5"
 
+[[Downloads]]
+deps = ["ArgTools", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+
 [[DrWatson]]
 deps = ["Dates", "FileIO", "LibGit2", "Pkg", "Random", "Requires", "UnPack"]
 git-tree-sha1 = "8a25e4984c3895ecc33891646362e4ff6baa4619"
@@ -281,10 +283,10 @@ uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 version = "6.2.0"
 
 [[GPUCompiler]]
-deps = ["DataStructures", "InteractiveUtils", "LLVM", "Libdl", "Scratch", "Serialization", "TimerOutputs", "UUIDs"]
-git-tree-sha1 = "c853c810b52a80f9aad79ab109207889e57f41ef"
+deps = ["DataStructures", "ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "Scratch", "Serialization", "TimerOutputs", "UUIDs"]
+git-tree-sha1 = "ef2839b063e158672583b9c09d2cf4876a8d3d55"
 uuid = "61eb1bfa-7361-4325-ad38-22787b887f55"
-version = "0.8.3"
+version = "0.10.0"
 
 [[Graphics]]
 deps = ["Colors", "LinearAlgebra", "NaNMath"]
@@ -450,7 +452,7 @@ version = "0.21.1"
 
 [[JUDI]]
 deps = ["ArgParse", "DSP", "Dierckx", "Distributed", "FFTW", "JOLI", "LinearAlgebra", "Pkg", "Printf", "PyCall", "Reexport", "SegyIO", "Test"]
-git-tree-sha1 = "16a467f456ffcae9d60be40e814c983b6b90bc4b"
+git-tree-sha1 = "48008ba9b6b2ae84caec1fe405576e31ba08da8a"
 repo-rev = "pyload"
 repo-url = "https://GitHub.com/slimgroup/JUDI.jl.git"
 uuid = "f3b833dc-6b2e-5b9c-b940-873ed6319979"
@@ -468,26 +470,24 @@ uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 version = "1.2.1"
 
 [[LazyArtifacts]]
-deps = ["Pkg"]
-git-tree-sha1 = "4bb5499a1fc437342ea9ab7e319ede5a457c0968"
+deps = ["Artifacts", "Pkg"]
 uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
-version = "1.3.0"
+
+[[LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
 
 [[LibCURL_jll]]
-deps = ["LibSSH2_jll", "Libdl", "MbedTLS_jll", "Pkg", "Zlib_jll", "nghttp2_jll"]
-git-tree-sha1 = "897d962c20031e6012bba7b3dcb7a667170dad17"
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
 uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
-version = "7.70.0+2"
 
 [[LibGit2]]
-deps = ["Printf"]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 
 [[LibSSH2_jll]]
-deps = ["Libdl", "MbedTLS_jll", "Pkg"]
-git-tree-sha1 = "717705533148132e5466f2924b9a3657b16158e8"
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
 uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
-version = "1.9.0+3"
 
 [[Libdl]]
 uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
@@ -540,10 +540,14 @@ deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
 [[MbedTLS_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "0eef589dd1c26a3ac9d753fe1a8bcad63f956fa6"
+deps = ["Artifacts", "Libdl"]
 uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
-version = "2.16.8+1"
+
+[[Memoize]]
+deps = ["MacroTools"]
+git-tree-sha1 = "2b1dfcba103de714d31c033b5dacc2e4a12c7caa"
+uuid = "c03570c3-d221-55d1-a50c-7939bbd78826"
+version = "0.4.4"
 
 [[Missings]]
 deps = ["DataAPI"]
@@ -566,6 +570,9 @@ git-tree-sha1 = "614e8d77264d20c1db83661daadfab38e8e4b77e"
 uuid = "e94cdb99-869f-56ef-bcf0-1ae2bcbe0389"
 version = "0.2.4"
 
+[[MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+
 [[NFFT]]
 deps = ["CUDA", "Distributed", "FFTW", "Graphics", "LinearAlgebra", "Random", "SparseArrays", "SpecialFunctions"]
 git-tree-sha1 = "44a74e77a99450c344edea4cae94d45a2f66ead8"
@@ -589,6 +596,9 @@ git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
 uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 version = "0.3.5"
 
+[[NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+
 [[Nullables]]
 git-tree-sha1 = "8f87854cc8f3685a60689d8edecaa29d2251979b"
 uuid = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd"
@@ -636,14 +646,14 @@ uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
 version = "1.1.0"
 
 [[Pkg]]
-deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
 [[Polynomials]]
 deps = ["Intervals", "LinearAlgebra", "OffsetArrays", "RecipesBase"]
-git-tree-sha1 = "1c6c5b0c3713738d6b987903c529d80622c37e07"
+git-tree-sha1 = "0b15f3597b01eb76764dd03c3c23d6679a3c32c8"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
-version = "1.2.0"
+version = "1.2.1"
 
 [[Primes]]
 deps = ["Test"]
@@ -668,7 +678,7 @@ uuid = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 version = "2.9.0"
 
 [[REPL]]
-deps = ["InteractiveUtils", "Markdown", "Sockets"]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 
 [[Random]]
@@ -769,9 +779,9 @@ version = "0.2.4"
 
 [[StaticArrays]]
 deps = ["LinearAlgebra", "Random", "Statistics"]
-git-tree-sha1 = "9da72ed50e94dbff92036da395275ed114e04d49"
+git-tree-sha1 = "2f01a51c23eed210ff4a1be102c4cc8236b66e5b"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.0.1"
+version = "1.1.0"
 
 [[Statistics]]
 deps = ["LinearAlgebra", "SparseArrays"]
@@ -783,8 +793,16 @@ git-tree-sha1 = "a83fa3021ac4c5a918582ec4721bc0cf70b495a9"
 uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 version = "0.33.4"
 
+[[TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[[Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+
 [[Test]]
-deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [[TextWrap]]
@@ -841,15 +859,13 @@ version = "0.5.3"
 
 [[XML2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
-git-tree-sha1 = "be0db24f70aae7e2b89f2f3092e93b8606d659a6"
+git-tree-sha1 = "afd2b541e8fd425cd3b7aa55932a257035ab4a70"
 uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
-version = "2.9.10+3"
+version = "2.9.11+0"
 
 [[Zlib_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
-git-tree-sha1 = "320228915c8debb12cb434c59057290f0834dbf6"
+deps = ["Libdl"]
 uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
-version = "1.2.11+18"
 
 [[Zstd_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -858,7 +874,9 @@ uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
 version = "1.4.8+0"
 
 [[nghttp2_jll]]
-deps = ["Libdl", "Pkg"]
-git-tree-sha1 = "8e2c44ab4d49ad9518f359ed8b62f83ba8beede4"
+deps = ["Artifacts", "Libdl"]
 uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
-version = "1.40.0+2"
+
+[[p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"

Project.toml

@@ -1,7 +1,7 @@
 name = "TimeProbeSeismic"
 uuid = "cd417b76-fc6b-4228-be91-94f5ea77264b"
 authors = ["Mathias Louboutin", "Felix J. Herrmann"]
-version = "0.0.1"
+version = "0.1.0"
 
 [deps]
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
@@ -21,5 +21,3 @@ Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 SlimOptim = "e4c7bc62-5b23-4522-a1b9-71c2be45f1df"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
-[compat]
-julia = "1.5.3"

scripts/splsrtm_seam.jl

@@ -14,10 +14,10 @@ vp, rho, n, d, o = deserialize(datadir("models", "bin.jld"))
 
 @warn "Untested script"
 
-m = vp.^(-2)
+m = 1/vp
 m0 = Float32.(imfilter(m, Kernel.gaussian(15)))
+m0 = m0.^2
 rho0 = Float32.(imfilter(rho[:, :], Kernel.gaussian(15)))
-dm = m0 - m
 
 n = size(m0)
 d = (12.5, 10.)
@@ -28,25 +28,18 @@ model0 = Model(n, d, o, m0, rho=rho0; nb=40)
 # Load data
 dobs, q = deserialize(datadir("data", "seam_obn_data.bin"))
 
-# Info structure for linear operators
-ntComp = get_computational_nt(q.geometry, dobs.geometry, model0)    # no. of computational time steps
-info = Info(prod(model0.n), dobs.nsrc, ntComp)
-
 #################################################################################################
 
-opt = Options(isic=false)
-
-# Setup operators
-F0 = judiModeling(info, model0, q.geometry, dobs.geometry; options=opt)
+opt = Options(isic=true, space_order=8)
 
 # Right-hand preconditioners (model topmute)
-idx_wb = find_water_bottom(reshape(dm, model0.n))
+idx_wb = find_water_bottom(vp .- minimum(vp))
 Tm = judiTopmute(model0.n, idx_wb, 10)  # Mute water column
 S = judiDepthScaling(model0)
 Mr = S*Tm
 
 # Curevelet transform
-C0 = joCurvelet2D(model0.n[1], 2*model0.n[2]; zero_finest = false, DDT = Float32, RDT = Float64)
+C0 = joCurvelet2D(model0.n[1], 2*model0.n[2]; zero_finest=false, DDT=Float32, RDT=Float32)
 function C_fwd(im, C, n)
     im = hcat(reshape(im, n), reshape(im, n)[:, end:-1:1])
     coeffs = C*vec(im)
@@ -61,70 +54,30 @@ end
 C = joLinearFunctionFwd_T(size(C0, 1), n[1]*n[2],
                           x -> C_fwd(x, C0, n),
                           b -> C_adj(b, C0, n),
-                          Float32,Float64, name="Cmirrorext")
-
-# Linearized Bregman parameters
-x = zeros(Float32, info.n)
-z = zeros(Float32, info.n)
-
-batchsize = 8
-niter = 5
-ps = 128
-fval = zeros(Float32, niter)
-tk = zeros(Float32, info.n)
+                          Float32,Float32, name="Cmirrorext")
 
-# Soft thresholding functions and Curvelet transform
-soft_thresholding(x::Array{Float64}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float64, lambda), 0.0)
-soft_thresholding(x::Array{Float32}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float32, lambda), 0f0)
-lambda = 0f0
-t = 1f-4
-
-# Main loop
+# Setup random shot selection
+batchsize = 4
 nsrc = 44
 ind_left = collect(1:nsrc÷2)
 ind_right = collect(nsrc÷2+1:nsrc)
 
-for j = 1:niter
+function breg_obj(x)
+    # Force log update
+    Base.flush(stdout)
     # Select batch and set up left-hand preconditioner
-    j == 1 ? batchs = batchsize : batchs = batchsize÷2
     if length(ind_left) == 0 || length(ind_right) == 0
         println("Made one pass through data, reseting source sampling")
         global ind_left = collect(1:nsrc÷2)
         global ind_right = collect(nsrc÷2+1:nsrc)
     end
-    i = vcat(sample_indices(ind_left, batchs), sample_indices(ind_right, batchs))
+    i = vcat(sample_indices(ind_left, batchsize), sample_indices(ind_right, batchsize))
+    batchsize == 4  && (global batchsize = 2)
     println("Iteration: $(j), imaging sources $(i)")
-    flush(Base.stdout)
-
-    phi, g = lsrtm_objective(model0, q[i], d_obs[i], Mr*x, ps; nlind=true, options=opt)
-
-    # Step size and update variable
-    g = adjoint(Mr)*g
-    fval[j] = phi
-
-    α = Float32(phi/norm(g)^2)
-    # sign
-    tk[:] .+= sign.(-g)
-
-    # Chatter correction
-    inds_z = findall(abs.(z) .> lambda)
-    stepg .= α
-    stepg[inds_z] .*= abs.(tk[inds_z])/j
-
-    @printf("At iteration %d function value is %2.2e and step length is %2.2e \n", j, fval[j], α)
-    flush(Base.stdout)
-    # Update variables and save snapshot
-    global z -= stepg[j, :] .* g
-
-    # Thresholding value
-    cz = C * z
-    (j-1)%10 == 0 && (global lambda = quantile(abs.(cz), .925))   # estimate thresholding parameter in 1st iteration
-    @printf("Lambda is %2.2e \n", lambda)
-
-    global x = adjoint(C)*soft_thresholding(cz, lambda)
-
-    savedict = @dict x z g lambda fval
-    wsave(datadir("lsrtm", "iter_$(j).bson"), savedict)
-end
 
+    f, g = lsrtm_objective(model0, q[i], dobs[i], Mr*x, ps;optins=opt, nlind=true)
+    return f, Mr'*g
+end
 
+opt = bregman_options(maxIter=20, verbose=2, quantile=.95, alpha=1, antichatter=true, spg=true)
+sol = bregman(breg_obj, 0f0.*vec(m0), opt, C)