<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 PDE
The acoustic wave equation for the square slowness m and a source q is given in 3D by 
Step1:  Time and space  discretization as a Taylor expansion.
The time discretization is define as a second order ( $ O (dt^2)) $) centered finite difference to get an explicit Euler scheme easy to solve by steping in time. 
$ \frac{d^2 u(x,t)}{dt^2} \simeq \frac{u(x,t+dt) - 2 u(x,t) + u(x,t-dt)}{dt^2} + O(dt^2) $
And we define the space discretization also as a Taylor serie, with oder chosen by the user. This can either be a direct expansion of the second derivative bulding the laplacian, or a combination of first oder space derivative. The second option can be a better choice in case you would want to extand the method to more complex wave equations involving first order derivatives in chain only.
$ \frac{d^2 u(x,t)}{dt^2} \simeq \frac{1}{dx^2} \sum_k \alpha_k (u(x+k dx,t)+u(x-k dx,t)) + O(dx^k) $
Step2:  Solve forward in time
The wave equation with absorbing boundary conditions writes
$ \eta \frac{d u(x,t)}{dt} + m \frac{d^2 u(x,t)}{dt^2} - \rho \nabla(\frac{1}{\rho} \text{grad}(u(x,t)) =q  $ 
and the adjont wave equation 
$ -\eta \frac{d u(x,t)}{dt} + m \frac{d^2 u(x,t)}{dt^2} - \rho \nabla(\frac{1}{\rho} \text{grad}(u(x,t)) =q  $ 
where $ \eta$  is a damping factor equal to zero inside the physical domain and decreasing inside the absorbing layer from the pysical domain to the border
And in order to simplify (kinda, I just like it that way) we will rewrite 
$ \rho \nabla(\frac{1}{\rho} \text{grad}(u(x,t)) = \nabla^2 u(x,t) + \rho \text{grad}(\frac{1}{\rho}) . \text{grad}(u(x,t)) $
Step3:  Define the discrete model
Step4:  Create functions for the PDE
The Gradient/Born are here so that everything is at the correct place, it is described later
Step5:  A Forward propagation example
Step6:  Adjoint test
In ordr to guaranty we have the gradient we need to make sure that the solution of the adjoint wave equation is indeed the true adjoint. Tod os so one should check that
$ <Ax,y> - <x,A^Ty> = 0$
where $A$ is the wave_equation,  $A^T$ is wave_equationA and $x,y$ are any random vectors in the range of each operator. This can however be expensive as this two vector would be of size $N * n_t$. To test our operator we will the relax this test by
$ <P_r A P_s^T x,y> - <x,P_SA^TP_r^Ty> = 0$
where $P_r , P_s^T$ are the source and recevier projection operator mapping the source and receiver locations and times onto the full domain. This allow to have only a random source of size $n_t$ at a random postion.
Step7:  Least square objective Gradient
We will consider here the least square objective, as this is the one in need of an adjoint. The test that will follow are however necessary for any objective and associated gradient in a optimization framework. The objective function can be written
$ min_m \Phi(m) 
Step8:  Adjoint test for the gradient
The adjoint of the FWI Gradient is the Born modelling operator, implementing a double propagation forward in time with a wavefield scaled by the model perturbation for the second propagation
$ J dm = - A^{-1}(\frac{d A^{-1}q}{dt^2}) dm $
Step9:  Jacobian test
The last part is to check that the operators are consistent with the problem. There is then two properties to be satisfied 
$ U(m + hdm) = U(m) +  \mathcal{O} (h) \
   U(m + h dm) = U(m) + h J[m]dm + \mathcal{O} (h^2) $
which are the linearization conditions for the objective. This is a bit slow to run here but here is the way to test it.
1 - Genrate data for the true model m
2 - Define a smooth initial model $m_0$ and comput the data $d_0$ for this model 
3 - You now have $U(m_0)$ 
4 - Define $ dm = m-m_0$ and $ h = {1,.1,.01,.001,...}$
5 - For each $h$ compute $U(m_0 + h dm)$ by generating data for $m_0 + h dm$ and compute $(J[m_0 + h dm]^T\delta |d) $ 
6 - Plot in Loglog the two lines of equation above
Step10:  Gradient test
The last part is to check that the operators are consistent with the problem. There is then two properties to be satisfied 
$ \Phi(m + hdm) = \Phi(m) +  \mathcal{O} (h) \
   \Phi(m + h dm) = \Phi(m) + h (J[m]^T\delta |d)dm + \mathcal{O} (h^2) $
which are the linearization conditions for the objective. This is a bit slow to run here but here is the way to test it.
1 - Genrate data for the true model m
2 - Define a smooth initial model $m_0$ and comput the data $d_0$ for this model 
3 - You now have $\Phi(m_0)$ 
4 - Define $ dm = m-m_0$ and $ h = {1,.1,.01,.001,...}$
5 - For each $h$ compute $\Phi(m_0 + h dm)$ by generating data for $m_0 + h dm$ and compute $(J[m_0 + h dm]^T\delta |d) $ 
6 - Plot in Loglog the two lines of equation above
<ASSISTANT_TASK:>
Python Code:
p=Function('p')
b=Function('b')
m,s,h,r = symbols('m s h r')
m=M(x,y)
q=Q(x,y,t)
d=D(x,y,t)
e=E(x,y)
r=rho(x,y)
Explanation: PDE
The acoustic wave equation for the square slowness m and a source q is given in 3D by :
\begin{cases} 
 &m \frac{d^2 u(x,t)}{dt^2} - \rho \nabla(\frac{1}{\rho} \text{grad}(u(x,t)) =q  \
 &u(.,0) = 0 \
 &\frac{d u(x,t)}{dt}|_{t=0} = 0 
\end{cases}
with the zero initial conditons to guaranty unicity of the solution
End of explanation
dtt=as_finite_diff(p(x,y,t).diff(t,t), [t-s,t, t+s])
dt=as_finite_diff(p(x,y,t).diff(t), [t-s, t+s])
# Spacial finite differences can easily be extended to higher order by increasing the list of sampling point in the next expression. 
# Be sure to keep this stencil symmetric and everything else in the notebook will follow.
dxx=as_finite_diff(p(x,y,t).diff(x,x), [x-h,x, x+h]) 
dyy=as_finite_diff(p(x,y,t).diff(y,y), [y-h,y, y+h])
dy=as_finite_diff(p(x,y,t).diff(y), [y-h, y+h])
dx=as_finite_diff(p(x,y,t).diff(x), [x-h, x+h])
dyr=as_finite_diff(b(x,y).diff(y), [y-h, y+h])
dxr=as_finite_diff(b(x,y).diff(x), [x-h, x+h])
dtt,dxx,dyy,dt,dx,dy
# In practice, compute rho. grad(rho) upfront (this is a constant) and just compute the following expression, it avoids 
# the recomputation of the gradient at every time step. int the following X and Y are the component of rho.grad(rho) in X and Y respectively
gradgrad=X*dx+Y*dy
expand(gradgrad)
lap=(dxx+dyy)
expand(lap + gradgrad)
Explanation: Time and space  discretization as a Taylor expansion.
The time discretization is define as a second order ( $ O (dt^2)) $) centered finite difference to get an explicit Euler scheme easy to solve by steping in time. 
$ \frac{d^2 u(x,t)}{dt^2} \simeq \frac{u(x,t+dt) - 2 u(x,t) + u(x,t-dt)}{dt^2} + O(dt^2) $
And we define the space discretization also as a Taylor serie, with oder chosen by the user. This can either be a direct expansion of the second derivative bulding the laplacian, or a combination of first oder space derivative. The second option can be a better choice in case you would want to extand the method to more complex wave equations involving first order derivatives in chain only.
$ \frac{d^2 u(x,t)}{dt^2} \simeq \frac{1}{dx^2} \sum_k \alpha_k (u(x+k dx,t)+u(x-k dx,t)) + O(dx^k) $
End of explanation
# Forward wave equation
wave_equation = m*dtt- (dxx+dyy)- r*(dxr*dx+dyr*dy) - q  + e*dt
stencil = solve(wave_equation,p(x,y,t+s))[0]
ts=lambdify((p(x,y,t-s),p(x-h,y,t), p(x,y,t), p(x+h,y,t),p(x,y-h,t), p(x,y+h,t),b(x-h,y), b(x,y), b(x+h,y),b(x,y-h), b(x,y+h), q , m,r, s, h,e),stencil,"numpy")
stencil
# Adjoint wave equation
wave_equationA =  m*dtt- (dxx+dyy) - r*(dxr*dx+dyr*dy)  - d  - e*dt
stencilA = solve(wave_equationA,p(x,y,t-s))[0]
tsA=lambdify((p(x,y,t+s),p(x-h,y,t), p(x,y,t), p(x+h,y,t),p(x,y-h,t), p(x,y+h,t),b(x-h,y), b(x,y), b(x+h,y),b(x,y-h), b(x,y+h), d, m,r, s, h,e),stencilA,"numpy")
stencilA
Explanation: Solve forward in time
The wave equation with absorbing boundary conditions writes
$ \eta \frac{d u(x,t)}{dt} + m \frac{d^2 u(x,t)}{dt^2} - \rho \nabla(\frac{1}{\rho} \text{grad}(u(x,t)) =q  $ 
and the adjont wave equation 
$ -\eta \frac{d u(x,t)}{dt} + m \frac{d^2 u(x,t)}{dt^2} - \rho \nabla(\frac{1}{\rho} \text{grad}(u(x,t)) =q  $ 
where $ \eta$  is a damping factor equal to zero inside the physical domain and decreasing inside the absorbing layer from the pysical domain to the border
And in order to simplify (kinda, I just like it that way) we will rewrite 
$ \rho \nabla(\frac{1}{\rho} \text{grad}(u(x,t)) = \nabla^2 u(x,t) + \rho \text{grad}(\frac{1}{\rho}) . \text{grad}(u(x,t)) $
End of explanation
import matplotlib.pyplot as plt
from matplotlib import animation
hstep=25 #space increment d  = minv/(10*f0);
tstep=2 #time increment dt < .5 * hstep /maxv;
tmin=0.0 #initial time
tmax=600 #simulate until
xmin=-500.0 - 10*hstep #left bound
xmax=500.0 + 10*hstep #right bound...assume packet never reaches boundary
ymin=-500.0 - 10*hstep #left bound
ymax=500.0 + 10*hstep #right bound...assume packet never reaches boundary
f0=.010
t0=1/.010
nbpml=10
nx = int((xmax-xmin)/hstep) + 1 #number of points on x grid
ny = int((ymax-ymin)/hstep) + 1 #number of points on x grid
nt = int((tmax-tmin)/tstep) + 2 #number of points on t grid
xsrc=-400
ysrc=0.0
xrec = nbpml+4
#set source as Ricker wavelet for f0 
def source(x,y,t):
    r = (np.pi*f0*(t-t0))
    val = (1-2.*r**2)*np.exp(-r**2)
    if abs(x-xsrc)<hstep/2 and abs(y-ysrc)<hstep/2:
        return val
    else:
        return 0.0
    
def dampx(x):
    dampcoeff=1.5*np.log(1.0/0.001)/(5.0*hstep);
    if x<nbpml:
        return dampcoeff*((nbpml-x)/nbpml)**2
    elif x>nx-nbpml-1:
        return dampcoeff*((x-nx+nbpml)/nbpml)**2
    else:
        return 0.0
def dampy(y):
    dampcoeff=1.5*np.log(1.0/0.001)/(5.0*hstep);
    if y<nbpml:
        return dampcoeff*((nbpml-y)/nbpml)**2
    elif y>ny-nbpml-1:
        return dampcoeff*((y-ny+nbpml)/nbpml)**2
    else:
        return 0.0
# Velocity and density models
def smooth10(vel,nx,ny):
    out=np.ones((nx,ny))
    out[:,:]=vel[:,:]
    for a in range(5,nx-6):
        out[a,:]=np.sum(vel[a-5:a+5,:], axis=0) /10
    return out
# True velocity
vel=np.ones((nx,ny)) + 2.0
rho=np.ones((nx,ny)) 
vel[floor(nx/2):nx,:]=4.5
rho[floor(nx/2):nx,:]=2.0
rho=rho**-1
mt=vel**-2
# Smooth velocity
v0=smooth10(vel,nx,ny)
m0=v0**-2
dm=m0-mt
Explanation: Define the discrete model
End of explanation
def Forward(nt,nx,ny,m,rho):
    u=np.zeros((nt,nx,ny))
    rec=np.zeros((nt,ny-2))
    for ti in range(0,nt):
        for a in range(1,nx-1):
            for b in range(1,ny-1):
                src = source(xmin+a*hstep,ymin+b*hstep,tstep*ti)
                damp=dampx(a)+dampy(b)
                if ti==0:
                    u[ti,a,b]=ts(0,0,0,0,0,0,rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
src,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                elif ti==1:
                    u[ti,a,b]=ts(0,u[ti-1,a-1,b],u[ti-1,a,b],u[ti-1,a+1,b],u[ti-1,a,b-1],u[ti-1,a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
src,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                else:
                    u[ti,a,b]=ts(u[ti-2,a,b],u[ti-1,a-1,b],u[ti-1,a,b],u[ti-1,a+1,b],u[ti-1,a,b-1],u[ti-1,a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
src,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                if a==xrec :
                    rec[ti,b-1]=u[ti,a,b] 
    return rec,u
def Adjoint(nt,nx,ny,m,rho,rec):
    v=np.zeros((nt,nx,ny))
    srca=np.zeros((nt))
    for ti in  range(nt-1, -1, -1):
        for a in range(1,nx-1):
            for b in range(1,ny-1):
                if a==xrec:
                    resid=rec[ti,b-1]
                else:
                    resid=0
                damp=dampx(a)+dampy(b)
                if ti==nt-1:
                    v[ti,a,b]=tsA(0,0,0,0,0,0,rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
resid,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                elif ti==nt-2:
                    v[ti,a,b]=tsA(0,v[ti+1,a-1,b],v[ti+1,a,b],v[ti+1,a+1,b],v[ti+1,a,b-1],v[ti+1,a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
resid,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                else:
                    v[ti,a,b]=tsA(v[ti+2,a,b],v[ti+1,a-1,b],v[ti+1,a,b],v[ti+1,a+1,b],v[ti+1,a,b-1],v[ti+1,a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
resid,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                if abs(xmin+a*hstep-xsrc)<hstep/2 and abs(ymin+b*hstep-ysrc)<hstep/2:
                    srca[ti]=v[ti,a,b]
    return srca,v
def Gradient(nt,nx,ny,m,rho,rec,u):
    v1=np.zeros((nx,ny))
    v2=np.zeros((nx,ny))
    v3=np.zeros((nx,ny))
    grad=np.zeros((nx,ny))
    for ti in range(nt-1,-1,-1):
        for a in range(1,nx-1):
            for b in range(1,ny-1):
                if a==xrec:
                    resid=rec[ti,b-1]
                else:
                    resid=0
                damp=dampx(a)+dampy(b)
                v3[a,b]=tsA(v1[a,b],v2[a-1,b],v2[a,b],v2[a+1,b],v2[a,b-1],v2[a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
resid,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                grad[a,b]=grad[a,b]-(v3[a,b]-2*v2[a,b]+v1[a,b])*(u[ti,a,b])
        v1,v2,v3=v2,v3,v1
    # No update inside the pml, only in the physical domain
    # grad[0:nbpml-1,:]=0
    # grad[nx-nbpml-1:nx-1,:]=0
    # grad[:,0:nbpml-1]=0
    # grad[:,ny-nbpml-1:ny-1]=0
    return tstep**-2*grad
def Born(nt,nx,ny,m,rho,dm):
    u1=np.zeros((nx,ny))
    U1=np.zeros((nx,ny))
    u2=np.zeros((nx,ny))
    U2=np.zeros((nx,ny))
    u3=np.zeros((nx,ny))
    U3=np.zeros((nx,ny))
    rec=np.zeros((nt,ny-2))
    src2=0
    for ti in range(0,nt):
        for a in range(1,nx-1):
            for b in range(1,ny-1):
                damp=dampx(a)+dampy(b)
                src = source(xmin+a*hstep,ymin+b*hstep,tstep*ti)
                u3[a,b]=ts(u1[a,b],u2[a-1,b],u2[a,b],u2[a+1,b],u2[a,b-1],u2[a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
src,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                src2 = -tstep**-2*(u3[a,b]-2*u2[a,b]+u1[a,b])*dm[a,b]
                U3[a,b]=ts(U1[a,b],U2[a-1,b],U2[a,b],U2[a+1,b],U2[a,b-1],U2[a,b+1],rho[a-1,b],rho[a,b],rho[a+1,b],rho[a,b-1],rho[a,b+1],
src2,m[a,b],rho[a,b]**-1,tstep,hstep,damp)
                if a==xrec :
                    rec[ti,b-1]=U3[a,b]
        u1,u2,u3=u2,u3,u1
        U1,U2,U3=U2,U3,U1
    return rec
Explanation: Create functions for the PDE
The Gradient/Born are here so that everything is at the correct place, it is described later
End of explanation
(rect,ut)=Forward(nt,nx,ny,mt,rho)
fig = plt.figure()
plts = []             # get ready to populate this list the Line artists to be plotted
plt.hold("off")
for i in range(nt):
    r = plt.imshow(ut[i,:,:])   # this is how you'd plot a single line...
    plts.append( [r] )  
ani = animation.ArtistAnimation(fig, plts, interval=50,  repeat = False)   # run the animation
#plt.show()
fig2 = plt.figure()
plt.hold("off")
shotrec = plt.imshow(rect,vmin=-10,vmax=10)   # this is how you'd plot a single line...
#plt.show()
Explanation: A Forward propagation example
End of explanation
(rec0,u0)=Forward(nt,nx,ny,m0,rho)
(srca,v)=Adjoint(nt,nx,ny,m0,rho,rec0)
plts = []             # get ready to populate this list the Line artists to be plotted
plt.hold("off")
for i in range(0,nt):
    r = plt.imshow(v[i,:,:])   # this is how you'd plot a single line...
    plts.append( [r] )  
ani = animation.ArtistAnimation(fig, plts, interval=50,  repeat = False)   # run the animation
#plt.show()
shotrec = plt.plot(srca)   # this is how you'd plot a single line...
plt.show()
# Actual adjoint test
term1=0
for ti in range(0,nt):
    term1=term1+srca[ti]*source(xsrc,ysrc,(ti)*tstep)
term2=LA.norm(rec0)**2
term1,term2,term1-term2,term1/term2
#if abs(term11/term21-1)<1e-9
#    print('Adjoint test passed')
Explanation: Adjoint test
In ordr to guaranty we have the gradient we need to make sure that the solution of the adjoint wave equation is indeed the true adjoint. Tod os so one should check that
$ <Ax,y> - <x,A^Ty> = 0$
where $A$ is the wave_equation,  $A^T$ is wave_equationA and $x,y$ are any random vectors in the range of each operator. This can however be expensive as this two vector would be of size $N * n_t$. To test our operator we will the relax this test by
$ <P_r A P_s^T x,y> - <x,P_SA^TP_r^Ty> = 0$
where $P_r , P_s^T$ are the source and recevier projection operator mapping the source and receiver locations and times onto the full domain. This allow to have only a random source of size $n_t$ at a random postion.
End of explanation
# Misfit
F0=.5*LA.norm(rec0-rect)**2
F0
Im1=Gradient(nt,nx,ny,m0,rho,rec0-rect,u0)
shotrec = plt.imshow(rec0-rect,vmin=-100,vmax=100)   # this is how you'd plot a single line...
#plt.show()
shotrec = plt.imshow(Im1,vmin=-100000,vmax=100000)   # this is how you'd plot a single line...
#plt.show()
Explanation: Least square objective Gradient
We will consider here the least square objective, as this is the one in need of an adjoint. The test that will follow are however necessary for any objective and associated gradient in a optimization framework. The objective function can be written
$ min_m \Phi(m) := \frac{1}{2} \| P_r A^{-1}(m) q - d\|_2^2$
And it's gradient becomes 
$ \nabla_m \Phi(m) = - (\frac{dA(m)u}{dm})^T v $
where v is the soltuion if the adjoint wave equation. For the simple acoustic case the  gradient can be rewritten as 
$ \nabla_m \Phi(m) = - \sum_{t=1}^{nt} \frac{d^2u(t)}{dt^2} v(t) $
End of explanation
Im2=Gradient(nt,nx,ny,m0,rho,rec0-rect,u0)
du1=Born(nt,nx,ny,m0,rho,dm)
term11=np.dot((rec0-rect).reshape(-1),du1.reshape(-1))
term21=np.dot(Im2.reshape(-1),dm.reshape(-1))
term11,term21,term11-term21,term11/term21
#if abs(term11/term21-1)<1e-9
#    print('Adjoint test passed')
shotrec = plt.imshow(Im2,vmin=-10000,vmax=10000)   # this is how you'd plot a single line...
#plt.show()
Explanation: Adjoint test for the gradient
The adjoint of the FWI Gradient is the Born modelling operator, implementing a double propagation forward in time with a wavefield scaled by the model perturbation for the second propagation
$ J dm = - A^{-1}(\frac{d A^{-1}q}{dt^2}) dm $
End of explanation
H=[1,0.1,0.01,.001,0.0001,0.00001,0.000001]
(D1,u0)=Forward(nt,nx,ny,m0,rho)
dub=Born(nt,nx,ny,m0,rho,dm)
error1=np.zeros((7))
error2=np.zeros((7))
for i in range(0,7):
    mloc=m0+H[i]*dm
    (d,u)=Forward(nt,nx,ny,mloc,rho)
    error1[i]  = LA.norm(d - D1,ord=1)
    error2[i]  = LA.norm(d - D1 - H[i]*dub,ord=1)
hh=np.zeros((7))
for i in range(0,7):
    hh[i]=H[i]*H[i]
shotrec = plt.loglog(H,error1,H,H)   # this is how you'd plot a single line...
plt.show()
shotrec = plt.loglog(H,error2,H,hh)   # this is howyou'd plot a single line...
plt.show()
Explanation: Jacobian test
The last part is to check that the operators are consistent with the problem. There is then two properties to be satisfied 
$ U(m + hdm) = U(m) +  \mathcal{O} (h) \
   U(m + h dm) = U(m) + h J[m]dm + \mathcal{O} (h^2) $
which are the linearization conditions for the objective. This is a bit slow to run here but here is the way to test it.
1 - Genrate data for the true model m
2 - Define a smooth initial model $m_0$ and comput the data $d_0$ for this model 
3 - You now have $U(m_0)$ 
4 - Define $ dm = m-m_0$ and $ h = {1,.1,.01,.001,...}$
5 - For each $h$ compute $U(m_0 + h dm)$ by generating data for $m_0 + h dm$ and compute $(J[m_0 + h dm]^T\delta |d) $ 
6 - Plot in Loglog the two lines of equation above
End of explanation
(DT,uT)=Forward(nt,nx,ny,mt,rho)
(D1,u0)=Forward(nt,nx,ny,m0,rho)
F0=.5*LA.norm(D1-DT)**2
g=Gradient(nt,nx,ny,m0,rho,D1-DT,u0)
G=np.dot(g.reshape(-1),dm.reshape(-1));
error21=np.zeros((7))
error22=np.zeros((7))
for i in range(0,7):
    mloc=m0+H[i]*dm
    (D,u)=Forward(nt,nx,ny,mloc,rho)
    error21[i]  = .5*LA.norm(D-DT)**2 -F0
    error22[i]  = .5*LA.norm(D-DT)**2 -F0 - H[i]*G
shotrec = plt.loglog(H,error21,H,H)   # this is how you'd plot a single line...
plt.show()
shotrec = plt.loglog(H,abs(error22),H,hh) # this is how you'd plot a single line...
plt.show()
error22
Explanation: Gradient test
The last part is to check that the operators are consistent with the problem. There is then two properties to be satisfied 
$ \Phi(m + hdm) = \Phi(m) +  \mathcal{O} (h) \
   \Phi(m + h dm) = \Phi(m) + h (J[m]^T\delta |d)dm + \mathcal{O} (h^2) $
which are the linearization conditions for the objective. This is a bit slow to run here but here is the way to test it.
1 - Genrate data for the true model m
2 - Define a smooth initial model $m_0$ and comput the data $d_0$ for this model 
3 - You now have $\Phi(m_0)$ 
4 - Define $ dm = m-m_0$ and $ h = {1,.1,.01,.001,...}$
5 - For each $h$ compute $\Phi(m_0 + h dm)$ by generating data for $m_0 + h dm$ and compute $(J[m_0 + h dm]^T\delta |d) $ 
6 - Plot in Loglog the two lines of equation above
End of explanation
<END_TASK>