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

.. only:: html

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

        :ref:`Go to the end <sphx_glr_download_examples_generated_scripts_synth_data_xray_tomography.py>`
        to download the full example code.

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

.. _sphx_glr_examples_generated_scripts_synth_data_xray_tomography.py:


Xray Tomography
===============

.. GENERATED FROM PYTHON SOURCE LINES 9-14

|Open In Colab|

.. |Open In Colab| image:: https://img.shields.io/badge/open%20in-Colab-b5e2fa?logo=googlecolab&style=flat-square&color=ffd670
   :target: https://colab.research.google.com/github/inlab-geo/cofi-examples/blob/main/examples/xray_tomography/xray_tomography.ipynb


.. GENERATED FROM PYTHON SOURCE LINES 17-24

If you are running this notebook locally, make sure you’ve followed
`steps
here <https://github.com/inlab-geo/cofi-examples#run-the-examples-with-cofi-locally>`__
to set up the environment. (This
`environment.yml <https://github.com/inlab-geo/cofi-examples/blob/main/envs/environment.yml>`__
file specifies a list of packages required to run the notebooks)


.. GENERATED FROM PYTHON SOURCE LINES 27-35

*Adapted from notebooks by Andrew Valentine & Malcolm Sambridge -
Research School of Earth Sciences, The Australian National University*

In this notebook, we look at an linear inverse problem based on Xray
Tomography. We will use ``cofi`` to run a linear system solver
(optionally with Tikhonov regularization and noise estimation) for this
problem.


.. GENERATED FROM PYTHON SOURCE LINES 38-44

0. Import modules
-----------------

The file ``xrt_tomography.py`` contains the forward code for this
problem.


.. GENERATED FROM PYTHON SOURCE LINES 44-55

.. code-block:: Python


    # -------------------------------------------------------- #
    #                                                          #
    #     Uncomment below to set up environment on "colab"     #
    #                                                          #
    # -------------------------------------------------------- #

    # !pip install -U cofi 
    # !git clone https://github.com/inlab-geo/cofi-examples.git
    # %cd cofi-examples/examples/xray_tomography








.. GENERATED FROM PYTHON SOURCE LINES 57-65

.. code-block:: Python


    import numpy as np
    import matplotlib.pyplot as plt

    from cofi import BaseProblem, InversionOptions, Inversion
    from cofi.utils import QuadraticReg
    import xrayTomography as xrt # import linear travel time forward model package








.. GENERATED FROM PYTHON SOURCE LINES 70-73

1. Define the problem
---------------------


.. GENERATED FROM PYTHON SOURCE LINES 73-82

.. code-block:: Python


    # display theory on the inference problem
    from IPython.display import display, Markdown

    with open("../../theory/geo_xray_tomography.md", "r") as f:
        content = f.read()

    display(Markdown(content))





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

 .. code-block:: none

    <IPython.core.display.Markdown object>




.. GENERATED FROM PYTHON SOURCE LINES 87-92

Firstly, we get some set up information for the problem. These include
the dataset and the Jacobian matrix. In the Xray Tomography example, the
Jacobian matrix is related to the lengths of paths within each grid.
Since the paths are fixed, the Jacobian matrix stays constant.


.. GENERATED FROM PYTHON SOURCE LINES 92-99

.. code-block:: Python


    #xrt = XrayTomography()
    # load data example
    loaded_dict = np.load("../../data/travel_time_tomography/linear_tomo_example.npz")
    linear_tomo_example = dict(loaded_dict)
    loaded_dict.close()








.. GENERATED FROM PYTHON SOURCE LINES 101-111

.. code-block:: Python


    # linear tomography problem set up
    paths = linear_tomo_example["_paths"]
    data = linear_tomo_example["_attns"]
    data_size = len(data)
    starting_model = linear_tomo_example["_start"]
    model_size,model_shape = starting_model.size,starting_model.shape
    good_model = linear_tomo_example["_true"]
    attns, jacobian = xrt.tracer(starting_model,paths)





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

 .. code-block:: none


    Evaluating paths:   0%|          | 0/10416 [00:00<?, ?it/s]
    Evaluating paths:   9%|▉         | 918/10416 [00:00<00:01, 9168.64it/s]
    Evaluating paths:  18%|█▊        | 1848/10416 [00:00<00:00, 9243.75it/s]
    Evaluating paths:  27%|██▋       | 2773/10416 [00:00<00:00, 9160.29it/s]
    Evaluating paths:  36%|███▌      | 3698/10416 [00:00<00:00, 9191.66it/s]
    Evaluating paths:  44%|████▍     | 4624/10416 [00:00<00:00, 9213.29it/s]
    Evaluating paths:  53%|█████▎    | 5546/10416 [00:00<00:00, 9186.60it/s]
    Evaluating paths:  62%|██████▏   | 6473/10416 [00:00<00:00, 9212.43it/s]
    Evaluating paths:  71%|███████   | 7398/10416 [00:00<00:00, 9222.99it/s]
    Evaluating paths:  80%|███████▉  | 8321/10416 [00:00<00:00, 9162.90it/s]
    Evaluating paths:  89%|████████▉ | 9255/10416 [00:01<00:00, 9214.58it/s]
    Evaluating paths:  98%|█████████▊| 10183/10416 [00:01<00:00, 9234.03it/s]
    Evaluating paths: 100%|██████████| 10416/10416 [00:01<00:00, 9215.58it/s]




.. GENERATED FROM PYTHON SOURCE LINES 113-118

.. code-block:: Python


    xrt_problem = BaseProblem()
    xrt_problem.set_data(data)
    xrt_problem.set_jacobian(jacobian)








.. GENERATED FROM PYTHON SOURCE LINES 123-126

We do some estimation on data noise and further perform a
regularization.


.. GENERATED FROM PYTHON SOURCE LINES 126-132

.. code-block:: Python


    sigma = 0.002
    #sigma = 0.1
    lamda = 50
    data_cov_inv = np.identity(data_size) * (1/sigma**2)








.. GENERATED FROM PYTHON SOURCE LINES 134-138

.. code-block:: Python


    xrt_problem.set_data_covariance_inv(data_cov_inv)
    xrt_problem.set_regularization(lamda * QuadraticReg(model_shape=(model_size,)))








.. GENERATED FROM PYTHON SOURCE LINES 143-145

Review what information is included in the ``BaseProblem`` object:


.. GENERATED FROM PYTHON SOURCE LINES 145-148

.. code-block:: Python


    xrt_problem.summary()





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

 .. code-block:: none

    =====================================================================
    Summary for inversion problem: BaseProblem
    =====================================================================
    Model shape: Unknown
    ---------------------------------------------------------------------
    List of functions/properties set by you:
    ['jacobian', 'regularization', 'data', 'data_covariance_inv']
    ---------------------------------------------------------------------
    List of functions/properties created based on what you have provided:
    ['jacobian_times_vector']
    ---------------------------------------------------------------------
    List of functions/properties that can be further set for the problem:
    ( not all of these may be relevant to your inversion workflow )
    ['objective', 'log_posterior', 'log_posterior_with_blobs', 'log_likelihood', 'log_prior', 'gradient', 'hessian', 'hessian_times_vector', 'residual', 'jacobian_times_vector', 'data_misfit', 'regularization_matrix', 'forward', 'data_covariance', 'initial_model', 'model_shape', 'blobs_dtype', 'bounds', 'constraints']




.. GENERATED FROM PYTHON SOURCE LINES 153-156

2. Define the inversion options
-------------------------------


.. GENERATED FROM PYTHON SOURCE LINES 156-160

.. code-block:: Python


    my_options = InversionOptions()
    my_options.set_tool("scipy.linalg.lstsq")








.. GENERATED FROM PYTHON SOURCE LINES 165-167

Review what’s been defined for the inversion we are about to run:


.. GENERATED FROM PYTHON SOURCE LINES 167-170

.. code-block:: Python


    my_options.summary()





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

 .. code-block:: none

    =============================
    Summary for inversion options
    =============================
    Solving method: None set
    Use `suggest_solving_methods()` to check available solving methods.
    -----------------------------
    Backend tool: `<class 'cofi.tools._scipy_lstsq.ScipyLstSq'>` - SciPy's wrapper function over LAPACK's linear least-squares solver, using 'gelsd', 'gelsy' (default), or 'gelss' as backend driver
    References: ['https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html', 'https://www.netlib.org/lapack/lug/node27.html']
    Use `suggest_tools()` to check available backend tools.
    -----------------------------
    Solver-specific parameters: None set
    Use `suggest_solver_params()` to check required/optional solver-specific parameters.




.. GENERATED FROM PYTHON SOURCE LINES 175-189

3. Start an inversion
---------------------

We can now solve the inverse problem using the Tikhonov-regularized form
of least-squares,

.. math:: \mathbf{m}=\left(\mathbf{A^TA}+\epsilon^2\sigma^2\mathbf{I}\right)^\mathbf{-1}\mathbf{A^Td}

where :math:`\sigma^2` is the variance of the expected noise on the
attenuation data.

For this dataset, we’ve taken :math:`\sigma = 0.002`\ s and chosen
:math:`\epsilon^2 = 50`.


.. GENERATED FROM PYTHON SOURCE LINES 189-194

.. code-block:: Python


    inv = Inversion(xrt_problem, my_options)
    inv_result = inv.run()
    inv_result.summary()





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

 .. code-block:: none

    ============================
    Summary for inversion result
    ============================
    SUCCESS
    ----------------------------
    model: [0.98494811 1.03000048 0.95776419 ... 0.94168322 1.03668701 1.00048943]
    sum_of_squared_residuals: []
    effective_rank: 2500
    singular_values: [9.30139732e+05 8.57631566e+05 8.57631566e+05 ... 1.14515274e+03
     8.80600410e+02 8.80600410e+02]
    model_covariance: [[ 1.17571588e-04 -8.57198189e-05 -1.62727362e-06 ...  1.56635037e-07
      -6.08653282e-08 -1.36217397e-07]
     [-8.57198189e-05  2.14596891e-04 -5.56362665e-05 ... -6.06195208e-07
       4.87748993e-07 -6.08653282e-08]
     [-1.62727362e-06 -5.56362665e-05  1.35540260e-04 ...  5.04358068e-07
      -6.06195208e-07  1.56635037e-07]
     ...
     [ 1.56635037e-07 -6.06195208e-07  5.04358068e-07 ...  1.35540260e-04
      -5.56362665e-05 -1.62727362e-06]
     [-6.08653282e-08  4.87748993e-07 -6.06195208e-07 ... -5.56362665e-05
       2.14596891e-04 -8.57198189e-05]
     [-1.36217397e-07 -6.08653282e-08  1.56635037e-07 ... -1.62727362e-06
      -8.57198189e-05  1.17571588e-04]]




.. GENERATED FROM PYTHON SOURCE LINES 199-205

4. Plotting
-----------

Below the two figures refers to the inferred model and true model
respectively.


.. GENERATED FROM PYTHON SOURCE LINES 205-209

.. code-block:: Python


    xrt.displayModel(inv_result.model.reshape(model_shape),clim=(1, 1.5),cmap=plt.cm.Blues); # inferred model
    xrt.displayModel(good_model,clim=(1, 1.5),cmap=plt.cm.Blues); # true model




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


    *

      .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_xray_tomography_001.png
         :alt: xray tomography
         :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_xray_tomography_001.png
         :class: sphx-glr-multi-img

    *

      .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_xray_tomography_002.png
         :alt: xray tomography
         :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_xray_tomography_002.png
         :class: sphx-glr-multi-img





.. GENERATED FROM PYTHON SOURCE LINES 214-231

5. Estimated uncertainties
--------------------------

We can now find the uncertainty on the recovered slowness parameters,
which describes how noise in the data propagate into the slowness
parameters with this data set. For the Tikhonov-regularised form of
least-squares, the model covariance matrix is a square matrix of size
:math:`M\times M`, where there are :math:`M` cells in the model.

.. math:: \mathbf{C_m}=\sigma^2\left(\mathbf{A^TA}+\epsilon^2\sigma^2\mathbf{I}\right)^\mathbf{-1}

.

This matrix was calculated as part of the solver routine above. The
square roots of the diagonal entries of this matrix are the
:math:`\sigma` errors in the slowness in each cell.


.. GENERATED FROM PYTHON SOURCE LINES 231-234

.. code-block:: Python


    Cm = inv_result.model_covariance








.. GENERATED FROM PYTHON SOURCE LINES 239-242

Lets plot the slowness uncertainties as a function of position across
the cellular model.


.. GENERATED FROM PYTHON SOURCE LINES 242-246

.. code-block:: Python


    slow_uncert = np.sqrt(np.diag(Cm)).reshape(model_shape)
    xrt.displayModel(slow_uncert,title='Slowness uncertainty',cmap=plt.cm.Greens)




.. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_xray_tomography_003.png
   :alt: Slowness uncertainty
   :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_xray_tomography_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 251-265

Uncertainty is uniformly low across the entire model and only
significant near the corners where there are few ray paths.

Similarly we can calculate uncertainty in velocity parameters using some
calculus.

.. math::  \Delta v = \left | \frac{\partial s}{\partial v}  \right | \Delta s 

and since :math:`s = 1/v` we get

.. math::  \Delta v = s^2\Delta s 

which gives the uncertainty image on velocity, which looks very similar.


.. GENERATED FROM PYTHON SOURCE LINES 265-268

.. code-block:: Python


    xrt.displayModel(slow_uncert * inv_result.model.reshape(model_shape),title='Velocity uncertainty',cmap=plt.cm.Blues);




.. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_xray_tomography_004.png
   :alt: Velocity uncertainty
   :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_xray_tomography_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 273-277

By clipping the colour range you can see an imprint of the true image,
indicating that high slowness/low velcoity areas have slightly higher
uncertainty.


.. GENERATED FROM PYTHON SOURCE LINES 280-293

--------------

Watermark
---------

.. raw:: html

   <!-- Feel free to add more modules in the watermark_list below, if more packages are used -->

.. raw:: html

   <!-- Otherwise please leave the below code cell unchanged -->


.. GENERATED FROM PYTHON SOURCE LINES 293-299

.. code-block:: Python


    watermark_list = ["cofi", "numpy", "scipy", "matplotlib"]
    for pkg in watermark_list:
        pkg_var = __import__(pkg)
        print(pkg, getattr(pkg_var, "__version__"))





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

 .. code-block:: none

    cofi 0.2.11
    numpy 2.3.5
    scipy 1.17.0
    matplotlib 3.10.8




.. GENERATED FROM PYTHON SOURCE LINES 305-305

sphinx_gallery_thumbnail_number = -1


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

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


.. _sphx_glr_download_examples_generated_scripts_synth_data_xray_tomography.py:

.. only:: html

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

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

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

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

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

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: xray_tomography.zip <xray_tomography.zip>`


.. only:: html

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

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