********************************
Part 2: The Adjoint State Method
********************************

In the previous section, you implemented the nonlinear forward operator

.. math::

    \mathcal{F}[m] = u.

In the `notes <http://math.mit.edu/icg/resources/notes325.pdf>`_, the
linearized version of this operator is derived and denoted by the operator
:math:`F`.  In this section, you will implement the adjoint of this operator,
:math:`F^*`, and see how the adjoint state method can be used to 1) produce
'images' from the imaging condition (Reverse Time Migration) and 2) compute
the gradient of the full waveform inversion nonlinear least squares objective
function.

The linear forward operator :math:`F` is addressed in full in the next
section.

Data Residuals
==============

The incident wavefield :math:`u_0(x,t)` is generated from the proposal
model :math:`c_0(x)`. The difference between the observed or measured
data :math:`d` and the simulated data :math:`d_0` due to :math:`u_0` is
called the residual :math:`r`, and

.. math::

   r(t) = d(t) - d_0(t) = d(t) - (\mathbf{S}u)(t) = d(t) - \mathbf{S}\mathcal{F}[m_0] (t).


.. topic:: Problem 2.1

    Generate the incident field :math:`u_0` from the proposal model :math:`m_0
    = \dfrac{1}{c_0^2}`. Compute and plot the residual trace. Be sure to use
    the same :math:`\Delta t` as you did for computing :math:`u(x,t)`,
    otherwise you will have to use interpolation to compute the residual.

Adjoint Wave Solver
===================

From the notes, we have that the computed model perturbation :math:`\delta m`,
which we will call the reverse time migration image :math:`I_\text{RTM}`, is
generated by applying the adjoint of the *linearized* forward model :math:`F`
to the residual,

.. math::

   I_\text{RTM} = F^*r,

or, when sampling is taken into account,

.. math::

   I_\text{RTM} = (\mathbf{S} F)^*r = F^{*} \mathbf{S}^{*} r = F^*r_{\text{ext}},

where :math:`r_\text{ext}` is the projection of the residual onto the full
wavefield.

Also from the notes,

.. math::

   F^*r_{\text{ext}} = -\int_0^T q(x,t) \partial_{tt}u_0(x,t) \textrm{d}t,

where :math:`q(x,t)` is the solution of

.. math::

   (\frac{1}{c_0(x)^2}\partial_{tt}-\partial_{xx})^{*}q_t(x,t) & = r_\text{ext}(x,t),  \\
   (\frac{1}{c_0(0)}\partial_t-\partial_x)^{*}q(0,t) & = 0, \\
   (\frac{1}{c_0(1)}\partial_t+\partial_x)^{*}q(1,t) & = 0, \\
   u(x,t) & = 0 \quad\text{for}\quad t \ge T,

where the adjoint operation :math:`*` is taken in time and collectively
:math:`L_0^{*}q(x,t) = r_\text{ext}(x,t)` (as opposed to
:math:`Lu(x,t)=f(x,t)` for the true wavefield). (Note that the argument
:math:`r_\text{ext}` can be replaced by any wavefield, :math:`F*` is defined
for any argument, not just the extended residual.)

.. topic:: Problem 2.2

    Extend your residual to the whole space using your ``point_source``
    function at the receiver location :math:`x_r`, to create a list of
    *adjoint source wavefields* . Using your ``leap_frog`` function, solve for
    the adjoint field :math:`q(x,t)`. Plot the space-time diagram of
    :math:`q(x,t)`.

    .. admonition:: Hint

        Be very careful, as the absorbing boundary conditions imply that
        :math:`L\ne L^*`. What does the adjoint operation mean in the time
        domain? How does this impact the time derivatives and your
        discretization of those derivatives?

Imaging Condition
=================

As mentioned before, the perturbation is

.. math::

   I_\text{RTM} = F^*\mathbf{S}^*r = \int_0^T\ q(x,t) \partial_{tt}u_0(x,t) \textrm{d}t,

this inner product is known as an *imaging condition* and is the core of a
gradient based optimization method. The image :math:`I_\text{RTM}` is said to
have been produced by *Reverse Time Migration*.

.. topic:: Problem 2.3

    Write a Python function ``imaging_condition(qs, u0s, config)`` that
    computes :math:`F^*r_\text{ext}`, given the incident wavefield
    :math:`u_0`.

    You may want to write a separate function to compute the second time
    derivative of the incident field. You can assume that :math:`u_0 = 0` for
    :math:`t > T`.

    .. admonition:: Hint

        Don't forget to account for the measure in your integrals.

    Plot the resulting image of the reflection. Compare that image to
    :math:`\delta C = C - C_0`. Think about how and why they are different?
    How might these results change if the source and receiver are on opposite
    sides of the reflector? What if there are multiple receivers?


    .. code:: python

        # Compute the image
        I_rtm = imaging_condition(qs, u0s, config)

        # Plot the comparison
        xs = np.arange(config['nx'])*config['dx']
        dC = C-C0

        plt.figure()
        plt.subplot(2, 1, 1)
        plt.plot(xs, dC, label=r'$\delta C$')
        plt.legend()
        plt.subplot(2, 1, 2)
        plt.plot(xs, I_rtm, label=r'$I_\text{RTM}$')
        plt.legend()

    .. image:: ../_static/exercises/p2_3.png

Visualizing the imaging condition as a space-time diagram is a very good way
to understand how the intersection of the incident and adjoint waves locates
the reflectors.

.. topic:: Problem 2.4

    Visualizing the imaging condition. Plot the space-time diagram for
    :math:`\partial_{tt}u_0(x,t)`. Either print the space-time diagrams of
    :math:`q(x,t)` and :math:`\partial_{tt}u_0(x,t)` and overlay them or use
    matplotlib to overlay the two images. See where the adjoint field and the
    incident field intersect to find the location of the reflection.  Try to
    understand why the intersections appear where they do.

Adjoint Operator
================

The adjoint operator is a linear operator and it will be useful to have a way
to apply it as a black box, much like we have a function ``forward_operator``
for applying the nonlinear forward operator as a black box.

.. topic:: Problem 2.5

    Write a function ``adjoint_operator(C0, d, config)`` which implements the
    adjoint operation you developed in Problems 2.2 and 2.3 and returns the
    reverse time migration image. Consider optionally returning the adjoint
    field.

    .. admonition:: Hint

        Remember that the argument ``d`` can be anything that 'looks like'
        data.  For example, the residual.
