.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "content/examples/20-published/plot_tomo_joint_with_volume.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_content_examples_20-published_plot_tomo_joint_with_volume.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_content_examples_20-published_plot_tomo_joint_with_volume.py:


Straight Ray with Volume Data Misfit Term
=========================================

Based on the SEG abstract Heagy, Cockett and Oldenburg, 2014.

Heagy, L. J., Cockett, A. R., & Oldenburg, D. W. (2014, August 5).
Parametrized Inversion Framework for Proppant Volume in a Hydraulically
Fractured Reservoir. SEG Technical Program Expanded Abstracts 2014.
Society of Exploration Geophysicists. doi:10.1190/segam2014-1639.1

This example is a simple joint inversion that consists of a

    - data misfit for the tomography problem
    - data misfit for the volume of the inclusions
      (uses the effective medium theory mapping)
    - model regularization

.. GENERATED FROM PYTHON SOURCE LINES 20-216



.. rst-class:: sphx-glr-horizontal


    *

      .. image-sg:: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_001.png
         :alt: Model Transform
         :srcset: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_001.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_002.png
         :alt: plot tomo joint with volume
         :srcset: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_002.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_003.png
         :alt: plot tomo joint with volume
         :srcset: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_003.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_004.png
         :alt: true, vol: 6.240e-01, recovered(no Volume term), vol: 2.652e-02 , recovered(with Volume term), vol: 3.795e-02 
         :srcset: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_004.png
         :class: sphx-glr-multi-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    True Volume: 0.6240000000000001

    Running inversion with SimPEG v0.22.0
    simpeg.InvProblem will set Regularization.reference_model to m0.
    simpeg.InvProblem will set Regularization.reference_model to m0.
    simpeg.InvProblem will set Regularization.reference_model to m0.
    simpeg.InvProblem will set Regularization.reference_model to m0.
    simpeg.InvProblem will set Regularization.reference_model to m0.

                        simpeg.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                        ***Done using the default solver Pardiso and no solver_opts.***
                    
    model has any nan: 0
    =============================== Projected GNCG ===============================
      #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   
    -----------------------------------------------------------------------------
    x0 has any nan: 0
       0  2.50e-01  8.48e+04  0.00e+00  8.48e+04    3.83e+01      0              
       1  2.50e-01  4.24e+03  9.97e-01  4.24e+03    3.67e+01      0              
    /home/vsts/work/1/s/simpeg/maps/_property_maps.py:1456: UserWarning:

    Maximum number of iterations reached

       2  2.50e-01  6.98e+02  8.59e-01  6.98e+02    3.24e+01      0              
       3  2.50e-01  6.67e+02  8.35e-01  6.67e+02    3.22e+01      4              
       4  2.50e-01  3.43e+02  7.28e-01  3.44e+02    2.57e+01      0              
       5  2.50e-01  8.63e+01  7.20e-01  8.64e+01    1.57e+01      0              
       6  2.50e-01  5.73e+01  7.23e-01  5.74e+01    1.51e+01      0              
       7  2.50e-01  4.99e+01  7.23e-01  5.00e+01    1.29e+01      0              
       8  2.50e-01  3.15e+01  7.26e-01  3.17e+01    1.10e+01      0              
       9  2.50e-01  2.31e+01  7.28e-01  2.32e+01    1.03e+01      0              
      10  2.50e-01  1.57e+01  7.32e-01  1.59e+01    8.89e+00      0              
      11  2.50e-01  1.26e+01  7.33e-01  1.28e+01    8.47e+00      0              
      12  2.50e-01  9.58e+00  7.36e-01  9.76e+00    7.66e+00      0              
      13  2.50e-01  6.81e+00  7.38e-01  6.99e+00    7.18e+00      0              
      14  2.50e-01  4.86e+00  7.38e-01  5.04e+00    6.23e+00      0              
      15  2.50e-01  3.60e+00  7.38e-01  3.78e+00    5.31e+00      0   Skip BFGS  
    ------------------------- STOP! -------------------------
    1 : |fc-fOld| = 1.2589e+00 <= tolF*(1+|f0|) = 8.4765e+03
    0 : |xc-x_last| = 4.5379e-01 <= tolX*(1+|x0|) = 1.0000e-01
    0 : |proj(x-g)-x|    = 5.3051e+00 <= tolG          = 1.0000e-01
    0 : |proj(x-g)-x|    = 5.3051e+00 <= 1e3*eps       = 1.0000e-02
    1 : maxIter   =      15    <= iter          =     15
    ------------------------- DONE! -------------------------

    Total recovered volume (no vol misfit term in inversion): 0.026515074507646754

    Running inversion with SimPEG v0.22.0

                        simpeg.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                        ***Done using the default solver Pardiso and no solver_opts.***
                    
    model has any nan: 0
    =============================== Projected GNCG ===============================
      #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment   
    -----------------------------------------------------------------------------
    x0 has any nan: 0
       0  2.50e-01  1.20e+05  0.00e+00  1.20e+05    1.56e+01      0              
       1  2.50e-01  7.79e+03  8.43e-01  7.79e+03    2.71e+01      0              
       2  2.50e-01  7.78e+03  8.15e-01  7.78e+03    2.65e+01      4              
       3  2.50e-01  7.77e+03  7.91e-01  7.78e+03    2.61e+01      4   Skip BFGS  
       4  2.50e-01  7.74e+03  7.57e-01  7.74e+03    2.55e+01      4              
       5  2.50e-01  7.64e+03  7.06e-01  7.64e+03    2.48e+01      3              
       6  2.50e-01  7.49e+03  6.73e-01  7.49e+03    2.39e+01      3              
    ------------------------------------------------------------------
    0 :    ft     = 7.4923e+03 <= alp*descent     = 7.4873e+03
    1 : maxIterLS =      10    <= iterLS          =     10
    ------------------------- End Linesearch -------------------------
    The linesearch got broken. Boo.

    Total volume (vol misfit term in inversion): 0.03795357151511156






|

.. code-block:: Python


    import numpy as np
    import scipy.sparse as sp
    import matplotlib.pyplot as plt

    from simpeg.seismic import straight_ray_tomography as tomo
    import discretize
    from simpeg import (
        maps,
        utils,
        regularization,
        optimization,
        inverse_problem,
        inversion,
        data_misfit,
        objective_function,
    )


    class Volume(objective_function.BaseObjectiveFunction):
        r"""
        A regularization on the volume integral of the model

        .. math::

            \phi_v = || \int_V m dV - \text{knownVolume} ||^2
        """

        def __init__(self, mesh, knownVolume=0.0, **kwargs):
            self.mesh = mesh
            self.knownVolume = knownVolume
            super().__init__(**kwargs)

        @property
        def knownVolume(self):
            """known volume"""
            return self._knownVolume

        @knownVolume.setter
        def knownVolume(self, value):
            self._knownVolume = utils.validate_float("knownVolume", value, min_val=0.0)

        def __call__(self, m):
            return (self.estVol(m) - self.knownVolume) ** 2

        def estVol(self, m):
            return np.inner(self.mesh.cell_volumes, m)

        def deriv(self, m):
            # return (self.mesh.cell_volumes * np.inner(self.mesh.cell_volumes, m))
            return (
                2
                * self.mesh.cell_volumes
                * (self.knownVolume - np.inner(self.mesh.cell_volumes, m))
            )  # factor of 2 from deriv of ||estVol - knownVol||^2

        def deriv2(self, m, v=None):
            if v is not None:
                return 2 * utils.mkvc(
                    self.mesh.cell_volumes * np.inner(self.mesh.cell_volumes, v)
                )
            else:
                # TODO: this is inefficent. It is a fully dense matrix
                return 2 * sp.csc_matrix(
                    np.outer(self.mesh.cell_volumes, self.mesh.cell_volumes)
                )


    def run(plotIt=True):
        nC = 40
        de = 1.0
        h = np.ones(nC) * de / nC
        M = discretize.TensorMesh([h, h])

        y = np.linspace(M.cell_centers_y[0], M.cell_centers_x[-1], int(np.floor(nC / 4)))
        rlocs = np.c_[0 * y + M.cell_centers_x[-1], y]
        rx = tomo.Rx(rlocs)

        source_list = [
            tomo.Src(location=np.r_[M.cell_centers_x[0], yi], receiver_list=[rx])
            for yi in y
        ]

        # phi model
        phi0 = 0
        phi1 = 0.65
        phitrue = utils.model_builder.create_block_in_wholespace(
            M.gridCC, [0.4, 0.6], [0.6, 0.4], [phi1, phi0]
        )

        knownVolume = np.sum(phitrue * M.cell_volumes)
        print("True Volume: {}".format(knownVolume))

        # Set up true conductivity model and plot the model transform
        sigma0 = np.exp(1)
        sigma1 = 1e4

        if plotIt:
            fig, ax = plt.subplots(1, 1)
            sigmaMapTest = maps.SelfConsistentEffectiveMedium(
                nP=1000, sigma0=sigma0, sigma1=sigma1, rel_tol=1e-1, maxIter=150
            )
            testphis = np.linspace(0.0, 1.0, 1000)

            sigetest = sigmaMapTest * testphis
            ax.semilogy(testphis, sigetest)
            ax.set_title("Model Transform")
            ax.set_xlabel(r"$\varphi$")
            ax.set_ylabel(r"$\sigma$")

        sigmaMap = maps.SelfConsistentEffectiveMedium(M, sigma0=sigma0, sigma1=sigma1)

        # scale the slowness so it is on a ~linear scale
        slownessMap = maps.LogMap(M) * sigmaMap

        # set up the problem and survey
        survey = tomo.Survey(source_list)
        problem = tomo.Simulation(M, survey=survey, slownessMap=slownessMap)

        if plotIt:
            _, ax = plt.subplots(1, 1)
            cb = plt.colorbar(M.plot_image(phitrue, ax=ax)[0], ax=ax)
            survey.plot(ax=ax)
            cb.set_label(r"$\varphi$")

        # get observed data
        data = problem.make_synthetic_data(phitrue, relative_error=0.03, add_noise=True)
        dpred = problem.dpred(np.zeros(M.nC))

        # objective function pieces
        reg = regularization.WeightedLeastSquares(M)
        dmis = data_misfit.L2DataMisfit(simulation=problem, data=data)
        dmisVol = Volume(mesh=M, knownVolume=knownVolume)
        beta = 0.25
        maxIter = 15

        # without the volume regularization
        opt = optimization.ProjectedGNCG(maxIter=maxIter, lower=0.0, upper=1.0)
        opt.remember("xc")
        invProb = inverse_problem.BaseInvProblem(dmis, reg, opt, beta=beta)
        inv = inversion.BaseInversion(invProb)

        mopt1 = inv.run(np.zeros(M.nC) + 1e-16)
        print(
            "\nTotal recovered volume (no vol misfit term in inversion): "
            "{}".format(dmisVol(mopt1))
        )

        # with the volume regularization
        vol_multiplier = 9e4
        reg2 = reg
        dmis2 = dmis + vol_multiplier * dmisVol
        opt2 = optimization.ProjectedGNCG(maxIter=maxIter, lower=0.0, upper=1.0)
        opt2.remember("xc")
        invProb2 = inverse_problem.BaseInvProblem(dmis2, reg2, opt2, beta=beta)
        inv2 = inversion.BaseInversion(invProb2)

        mopt2 = inv2.run(np.zeros(M.nC) + 1e-16)
        print("\nTotal volume (vol misfit term in inversion): {}".format(dmisVol(mopt2)))

        # plot results

        if plotIt:
            fig, ax = plt.subplots(1, 1)
            ax.plot(data.dobs)
            ax.plot(dpred)
            ax.plot(problem.dpred(mopt1), "o")
            ax.plot(problem.dpred(mopt2), "s")
            ax.legend(["dobs", "dpred0", "dpred w/o Vol", "dpred with Vol"])

            fig, ax = plt.subplots(1, 3, figsize=(16, 4))
            im0 = M.plot_image(phitrue, ax=ax[0])[0]
            im1 = M.plot_image(mopt1, ax=ax[1])[0]
            im2 = M.plot_image(mopt2, ax=ax[2])[0]

            for im in [im0, im1, im2]:
                im.set_clim([0.0, phi1])

            plt.colorbar(im0, ax=ax[0])
            plt.colorbar(im1, ax=ax[1])
            plt.colorbar(im2, ax=ax[2])

            ax[0].set_title("true, vol: {:1.3e}".format(knownVolume))
            ax[1].set_title(
                "recovered(no Volume term), vol: {:1.3e} ".format(dmisVol(mopt1))
            )
            ax[2].set_title(
                "recovered(with Volume term), vol: {:1.3e} ".format(dmisVol(mopt2))
            )

            plt.tight_layout()


    if __name__ == "__main__":
        run()
        plt.show()


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 35.356 seconds)

**Estimated memory usage:**  9 MB


.. _sphx_glr_download_content_examples_20-published_plot_tomo_joint_with_volume.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_tomo_joint_with_volume.ipynb <plot_tomo_joint_with_volume.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_tomo_joint_with_volume.py <plot_tomo_joint_with_volume.py>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_