velocoti gradient

.github/workflows/CI.yml

@@ -35,7 +35,7 @@ jobs:
 
     steps:
       - name: Checkout PhotoAcoustic
-        uses: actions/checkout@v2
+        uses: actions/checkout@v3
 
       - name: Setup julia
         uses: julia-actions/setup-julia@v1

.github/workflows/docs.yml

@@ -10,7 +10,7 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
 
       - uses: julia-actions/setup-julia@latest
 

Project.toml

@@ -1,7 +1,7 @@
 name = "PhotoAcoustic"
 uuid = "86b14aa7-fcb7-4836-b4c7-056f45a9c77b"
 authors = ["rafaelorozco <[email protected]>", "Mathias Louboutin <[email protected]"]
-version = "0.3.0"
+version = "0.4.0"
 
 [deps]
 JUDI = "f3b833dc-6b2e-5b9c-b940-873ed6319979"
@@ -10,7 +10,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-JUDI = "^3.1.3"
+JUDI = "^3.1.10"
 Reexport = "1"
 julia = "1.6"
 

docs/src/derivations.md

@@ -86,7 +86,15 @@ We derive the sensitivity with repect to the square slowness in a similar fation
 
 
 ```math
-    d_{m}F = - \int_{0}^{T} \lambda^{\top} \ddot u
+    d_{m}F = - \int_{0}^{T} \lambda^{\top} \ddot u dt
 ```
 
-where $\lambda$ is the solution of the same adjoint equation derived for the sensitivity with respect to the photoacoustic source.
+where $\lambda$ is the solution of the same adjoint equation derived for the sensitivity with respect to the photoacoustic source. Note that unlike the source excitation case (i.e point source in [JUDI]) this is not strictly equivalent to $\ddot \lambda^{\top} u$ due to the intial state. Because JUDI computes $- \int_{0}^{T} \ddot \lambda^{\top} u dt$ the following correction is needed:
+
+
+```math
+\begin{aligned}
+    d_{m}F &= - \int_{0}^{T} \lambda^{\top} \ddot u dt =  - \int_{0}^{T}  \ddot \lambda^{\top} u dt  + \lambda(0)  \dot u (0) - \dot \lambda(0) u(0) \\
+    &=  - \int_{0}^{T}  \ddot \lambda^{\top} u dt - \dot \lambda(0) p_0
+\end{aligned}
+```

src/implementation.py

@@ -81,20 +81,19 @@ def adjointbornis(model, y, rcv_coords, init_dist, **kwargs):
     """
     nt = y.shape[0]
     born_fwd = kwargs.get('born_fwd', False)
-    rec, u, _ = op_fwd_JIS[born_fwd](model, rcv_coords, init_dist, nt, **kwargs)
+    rec, u, _ = op_fwd_JIS[born_fwd](model, rcv_coords, init_dist, nt,
+                                     save=True, **kwargs)
 
-    kwargs['return_op'] = True
     kwargs.pop('freq_list', None)
-    op, g, kwg = gradient(model, -y, rcv_coords, u, **kwargs)
+    kwargs['return_op'] = True
+    op, g, kwg = gradient(model, y, rcv_coords, u, save=True, **kwargs)
     op(**kwg)
 
-    # Need the extra v[0].dt * m
-    # v = kwg['v']
-    # Correct for default scaling in injection
-    # mrm = model.m * model.irho
-    # from IPython import embed; embed()
-    # op = Operator(Eq(g, g / mrm))
-    # op(dt=model.critical_dt, time_m=0, time_M=0)
+    # Need the intergation by part correction since we compute the gradient on
+    # u * v.dt (see Documentation)
+    v = kwg['v']
+    op = Operator(Eq(g, g - v.dt*u))
+    op(dt=model.critical_dt, time_m=0, time_M=0)
 
     return g.data
 

src/judiInitialState.jl

@@ -30,6 +30,7 @@ Arguments
 judiInitialState(source::Vector{Array{T, N}}) where {T<:Number, N} = judiInitialState{T}(length(source), source)
 judiInitialState(source::Array{T, N}) where {T<:Number, N} = judiInitialState([source])
 judiInitialState(source::Array{T, N}, nsrc::Integer) where {T<:Number, N} = judiInitialState([source for s=1:nsrc])
+judiInitialState(x::judiInitialState) = x
 
 ############################################################
 # JOLI conversion

src/judiPhoto.jl

@@ -66,7 +66,7 @@ function _forward_prop(J::judiPhoto{T, O}, q::AbstractArray{T}, op::PyObject; dm
     nt = length(0:dtComp:recGeometry.t[1])
 
     # Set up coordinates
-    rec_coords = setup_grid(recGeometry, J.F.model.n)    # shifts rec coordinates by origin
+    rec_coords = setup_grid(recGeometry, J.F.model.n)
 
     # Devito interface
     dsim = wrapcall_data(op, modelPy, rec_coords, init_dist, nt, space_order=J.F.options.space_order)
@@ -87,7 +87,7 @@ function _reverse_propagate(J::judiPhoto{T, O}, q::AbstractArray{T}, op::PyObjec
     dtComp  = convert(Float32, modelPy."critical_dt")
 
     # Set up coordinates
-    rec_coords = setup_grid(recGeometry, J.F.model.n)    # shifts rec coordinates by origin
+    rec_coords = setup_grid(recGeometry, J.F.model.n)
 
     # Extrapolate input data to computational grid     
     qIn = time_resample(srcData, recGeometry, dtComp)[1]

test/runtests.jl

@@ -2,7 +2,6 @@ using PhotoAcoustic
 
 using LinearAlgebra, Test, Printf
 
-
 ######## Copy paste from JUDI, should reorganize JUDI so can just be imported
 test_adjoint(f::Bool, j::Bool, last::Bool) = (test_adjoint(f, last), test_adjoint(j, last))
 test_adjoint(adj::Bool, last::Bool) = (adj || last) ? (@test adj) : (@test_skip adj)
@@ -58,9 +57,9 @@ function grad_test(misfit, x0, dx, g; maxiter=6, h0=5f-2, data=false, stol=1f-1)
 
     rate1 = err1[1:end-1]./err1[2:end]
     rate2 = err2[1:end-1]./err2[2:end]
-    # @test isapprox(mean(rate1), 1.25f0; atol=stol)
-    # @test isapprox(mean(rate2), 1.5625f0; atol=stol)
+    @test isapprox(mean(rate1), 1.25f0; atol=stol)
+    @test isapprox(mean(rate2), 1.5625f0; atol=stol)
 end
 
 
-# include("test_adjoints.jl")
+include("test_sensitivities.jl")

test/test_adjoints.jl → test/test_sensitivities.jl

@@ -2,8 +2,6 @@ using PhotoAcoustic
 
 using LinearAlgebra, Test, Printf
 
-include("./runtests.jl")
-
 # Set up model structure
 n = (80, 80)   # (x,y,z) or (x,z)
 d = (0.08f0, 0.08f0)
@@ -25,9 +23,9 @@ model0 = Model(n, d, o, m0;)
 
 # Set up receiver geometry
 nxrec = 64
-xrec = range(0, stop=d[1]*(n[1]-1), length=nxrec)
+xrec = range(0.08f0, stop=d[1]*(n[1]-2), length=nxrec)
 yrec = [0f0]
-zrec = range(0, stop=0, length=nxrec)
+zrec = range(0.08f0, stop=0.08f0, length=nxrec)
 
 # receiver sampling and recording time
 time = 5f0 #[microsec] 
@@ -45,33 +43,28 @@ tol = 5f-4
 maxtry = 3
 
 # Options
-opt = Options(sum_padding=true, dt_comp=dt)
+opt = Options(sum_padding=true)
+w = zeros(Float32, model.n...)
+w[35:45, 35:45] .= randn(Float32, 11, 11)
+w = judiInitialState(w)
+
+F = judiPhoto(model, recGeometry; options=opt)
+dobs = F*w
 
 @testset "Jacobian test" begin
-    w = zeros(Float32, model.n...)
-    w[35:45, 35:45] .= randn(Float32, 11, 11)
-    w = judiInitialState(w)
+
     # Setup operators
-    F = judiPhoto(model, recGeometry; options=opt)
     F0 = judiPhoto(model0, recGeometry; options=opt)
     J = judiJacobian(F0, w)
 
     # Linear modeling
-    dobs = F*w
     dD = J*vec(dm)
 
     # Gradient test
     grad_test(x-> F(;m=x)*w, m0, dm, dD; data=true)
 end
 
 @testset "FWI gradient test" begin
-    w = zeros(Float32, model.n...)
-    w[35:45, 35:45] .= randn(Float32, 11, 11)
-    w = judiInitialState(w)
-    # Setup operators
-    F = judiPhoto(model, recGeometry; options=opt)
-    dobs = F*w
-
     # Background operators
     F0 = judiPhoto(model0, recGeometry; options=opt)
     J = judiJacobian(F0, w)
@@ -83,19 +76,15 @@ end
 end
 
 @testset "Photoacoustic  adjoint test with constant background" begin
-    F = judiModeling(model; options=opt)
+    Ainv = judiModeling(model; options=opt)
     Pr = judiProjection(recGeometry)
     I = judiInitialStateProjection(model)
 
 	# Setup operators
-	A = judiPhoto(F, recGeometry)
-    A2 = Pr*F*I'
+	A = judiPhoto(Ainv, recGeometry)
+    A2 = Pr*Ainv*I'
     A3 = judiPhoto(model, recGeometry; options=opt)
 
-    w = zeros(Float32, model.n...)
-    w[35:45, 35:45] .= randn(Float32, 11, 11)
-    w = judiInitialState(w)
-
     # Nonlinear modeling
     y = A*w
     y2 = A2*w