lsrtm_objective does not produce the same gradient as is produced by Jacobians

[ziyiyin97]

I compared ways to get gradients in LS-RTM

1. Use lsrtm_objective.
2. Do linear algebras by judiJacobian operator.

and they produce different results. However, there is no problem when x is all 0. The script is shown below.

Also a small thing: J[1] works by J[[1]] produces error at https://github.com/slimgroup/JUDI.jl/blob/dda83aa9971002241fe9e0fed6285f5f21cb7547/src/TimeModeling/Types/GeometryStructure.jl#L294 because of unsupported keyword arguments.

using PyPlot
using Random, Images, JLD2, LinearAlgebra
using JOLI, Statistics
using JUDI
Random.seed!(1234)

n = (101,101)
d = (12.5f0,12.5f0)
o = (0f0, 0f0)
v = 1.5*ones(Float32,n)
v[:,19:end] .= 2.5f0
v[:,40:end] .= 3f0
v[:,60:end] .= 3.5f0
m = 1f0./v.^2f0
m0 = convert(Array{Float32,2},imfilter(m, Kernel.gaussian(4)))

model0 = Model(n,d,o,m0;nb=80)

nsrc = 1
nrec = n[1]

dtS = dtR = 1f0
timeS = timeR = 2000f0

xsrc = convertToCell(range((n[1]-1)*d[1]/2,stop=(n[1]-1)*d[1]/2,length=nsrc))
ysrc = convertToCell(range(0f0,stop=0f0,length=nsrc))
zsrc = convertToCell(range(17*d[2]-2f0,stop=17*d[2]-2f0,length=nsrc))

srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)
wavelet = ricker_wavelet(timeS,dtS,0.012f0)

q = judiVector(srcGeometry, wavelet)

xrec = range(d[1],stop=(n[1]-1)*d[1],length=nrec)
yrec = 0f0
zrec = range(10f0,stop=10f0,length=nrec)

recGeometry = Geometry(xrec,yrec,zrec; dt=dtR, t=timeR, nsrc=nsrc)

ntComp = get_computational_nt(srcGeometry, recGeometry, model0)
info = Info(prod(n), nsrc, ntComp)

opt = Options(isic=true)

F0 = judiProjection(info, recGeometry)*judiModeling(info, model0; options=opt)*adjoint(judiProjection(info, srcGeometry))
J = judiJacobian(F0,q)
dlin = J*vec(m-m0)

x = J'*dlin

phi1, g1 = lsrtm_objective(model0, q, dlin, x; nlind=false,options=opt)
g1 = vec(g1)
d_res = J*x-dlin
g = J'*d_res
g = vec(g.data)
phi = 0.5f0 * norm(d_res)^2

println("misfit in g = ", norm(g1-g)/norm(g1+g)) ## problem here
println("misfit in phi = ", norm(phi1-phi)/norm(phi1+phi)) ## problem here

x1 = zeros(Float32, prod(n))

phi1, g1 = lsrtm_objective(model0, q, dlin, x1; nlind=false,options=opt)
g1 = vec(g1)
d_res = J*x1-dlin
g = J'*d_res
g = vec(g.data)
phi = 0.5f0 * norm(d_res)^2

println("misfit in g = ", norm(g1-g)/norm(g1+g))    ## no problem here
println("misfit in phi = ", norm(phi1-phi)/norm(phi1+phi))  ## no problem here

[ziyiyin97]

This error only appears when isic=true

[ziyiyin97]

Problem at https://github.com/slimgroup/JUDI.jl/blob/64c2a7d0d988eca4ff541f6b09b85e207c4b26f8/src/pysource/interface.py#L585 Will be updated in PR