Probing judiJacobian only does single source

[ziyiyin97]

MFE is attached. adjoint(judiJacobian(F0,q,ps,dD))*dD seems only to do RTM probing on the first source (while F0,q include multiple sources). lsrtm_objective works well!

using Pkg
Pkg.activate("TimeProbeSeismic")
Pkg.add(url="https://github.com/slimgroup/JUDI.jl.git",rev="pyload")
Pkg.update()
Pkg.instantiate()
using DrWatson
@quickactivate :TimeProbeSeismic
using JLD2, FFTW, DSP

# Set up model structure
n = (1200, 200)   # (x,y,z) or (x,z)
d = (10., 10.)
o = (0., 0.)

# Velocity [km/s]
v = ones(Float32,n) .+ 0.5f0
v[:,30:end] .= 3.5f0

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2
m0 = convert(Array{Float32,2},imfilter(m, Kernel.gaussian(3)))
dm = vec(m - m0)

# Setup info and model structure
nsrc = 2	# number of sources
model = Model(n, d, o, m)
model0 = Model(n, d, o, m0)

# Set up receiver geometry
nxrec = 480
xrec = range(3000f0, stop=9000f0, length=nxrec)
yrec = 0f0
zrec = range(50f0, stop=50f0, length=nxrec)

# receiver sampling and recording time
timeR = 1000f0   # receiver recording time [ms]
dtR = 2f0    # receiver sampling interval [ms]

# Set up receiver structure
recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

# Set up source geometry (cell array with source locations for each shot)
xsrc = convertToCell(range(5000f0, stop=7000f0, length=nsrc))
ysrc = convertToCell(range(0f0, stop=0f0, length=nsrc))
zsrc = convertToCell(range(200f0, stop=200f0, length=nsrc))

# source sampling and number of time steps
timeS = 1000f0  # ms
dtS = 2f0   # ms

# Set up source structure
srcGeometry = Geometry(xsrc, ysrc, zsrc; dt=dtS, t=timeS)

# setup wavelet
f0 = 0.015f0     # kHz
wavelet = ricker_wavelet(timeS, dtS, f0)
q = judiVector(srcGeometry, wavelet)

# Set up info structure for linear operators
ntComp = get_computational_nt(srcGeometry, recGeometry, model)
info = Info(prod(n), nsrc, ntComp)

###################################################################################################

# Write shots as segy files to disk
opt = Options(isic=true)

# Setup operators
Pr = judiProjection(info, recGeometry)
F0 = judiModeling(info, model0; options=opt)
Ps = judiProjection(info, srcGeometry)
F0 = Pr*F0*adjoint(Ps)
J = judiJacobian(F0, q)

# Linearized modeling
dD = J*dm
# Adjoint jacobian
println("conventional RTM")
@time rtm = adjoint(J)*dD
# probing RTM

println("32 probing vectors")
@time rtm_32 = adjoint(judiJacobian(F0,q,32,dD))*dD

figure();imshow(rtm.data',aspect=3,cmap="Greys",vmin=-0.1*norm(rtm.data,Inf),vmax=0.1*norm(rtm.data,Inf));title("conventional RTM")
figure();imshow(rtm_32.data',aspect=3,cmap="Greys",vmin=-0.1*norm(rtm.data,Inf),vmax=0.1*norm(rtm.data,Inf));title("RTM w/ 32 probing vectors")

f, g = JUDI.lsrtm_objective(model0, q, dD, 0f0 .* dm; nlind=false, options=opt)
f_32, g_32 = JUDI.lsrtm_objective(model0, q, dD, 0f0 .* dm, 32; nlind=false, options=opt)

figure();imshow(-g',aspect=3,cmap="Greys",vmin=-0.1*norm(rtm.data,Inf),vmax=0.1*norm(rtm.data,Inf));title("conventional RTM (by lsrtm_objective)")
figure();imshow(-g_32',aspect=3,cmap="Greys",vmin=-0.1*norm(rtm.data,Inf),vmax=0.1*norm(rtm.data,Inf));title("RTM w/ 32 probing vectors (by lsrtm_objective)")

[mloubout]

srcnum missing in adjbornop is one line fix